text:
# Extinction calculations of multi-sphere polycrystalline graphitic clusters

## 1 Introduction

The shape of interstellar and circumstellar grains is still an outstanding issue. For many years, the complexity of the electromagnetic scattering problem to solve has limited the shapes studied to spheres, infinite cylinders and spheroids. However, the shape of many interstellar grains are expected to be non-spherical and maybe even highly irregular. One way to deal with irregular particles and so with clusters of dust grains is to assume that they consist of touching spheres. With such an assumption it is possible to construct many distinct different morphologies which can then be compared with observations.

The problem of evaluating the extinction efficiency ($`Q_{\mathrm{ext}}`$) is that of solving Maxwell’s equations with appropriate boundary conditions at the cluster surface. A solution was formulated by Lorenz (lorenz (1890)) and Mie (mie (1908)) for a homogeneous single sphere and the complete formalism is therefore often referred to as Lorenz-Mie theory. A complementary solution based on expansion of scalar potentials was given by Debye (debye (1909)). A detailed description of this exact electromagnetic solution can be found in Bohren & Huffman (bohren+huffman (1983)). For a comprehensive review on the optics of cosmic dust see Voshchinnikov (nicolai02 (2002)).

An updated overview of exact theories and numerical techniques for computing the scattered electromagnetic field by clusters of particles is given in (Mishchenko et al. 2000a ; 2000b ; Fuller & Mackowski Fuller (2000); Ciric & Coorey Ciric (2000); Draine draine00 (2000); Voshchinnikov nicolai02 (2002) and references therein). All of these methods are based on solving Maxwell’s equations. For clusters of spheres embedded in non-absorbing media one such method is based on the Gérardy and Ausloos theory (Gérardy & Ausloos GA4 (1980); GA (1982); GA2 (1983); GA3 (1984)) which recently has been extended to also treat clusters embedded in absorbing media (Lebedev & Stenzel lebedev1 (1999); Lebedev et al. lebedev2 (1999)). Another method that is often used in practice is the discrete dipole approximation (DDA); it has been used in a wide range of scattering problems concerning clusters of particles including the extinction of interstellar dust grains (e.g. Vaidya et al. vaidya (2001); Bazell & Dwek bazell (1990)).

The Galactic interstellar extinction curve displays a “2175 Å peak”, which has presented an astrophysical puzzle since its discovery by Stecher (stecher (1965)). Fitzpatrick & Massa (fitzpatrick+massa86 (1986)) studied the interstellar extinction in the direction of 45 reddened stars, and found that it displays a peak whose central wavelength $`\lambda _0`$ is remarkably constant ($`\lambda _0=2174.4\pm 17`$ Å), even though its full width at half maximum (FWHM) varies considerably from $`360`$ to $`600`$ Å; they also found no apparent correlation between the small variation in $`\lambda _0`$ and the large variation in the FWHM.

Graphite has long been considered a very promising though controversial candidate for explaining the 2175 Å peak (e.g. Stecher & Donn stecher+donn (1965); Fitzpatrick & Massa fitzpatrick+massa (1988); Mathis mathis (1994); Voshchinnikov nicolai90 (1990); Sorell sorell (1990); Draine & Malhotra draine+malhotra (1993); Rouleau et al. rouleau (1997); Will & Aannestad will+aannestad (1999)). It is considered a promising candidate because small, 0.015 $`\mu `$m, graphitic spheres display an extinction peak at about the right wavelength with a FWHM in accordance with the observations (Gilra gilra (1972)), and because its abundance do not seem to contradict the cosmic abundance constraints (Snow & Witt snow+witt (1995)). It is considered a controversial candidate, however, because increasing the size of a small graphitic particle simultaneously increases the central wavelength and FWHM of its extinction peak. The peak will shift to longer wavelengths when the particle size is increased, when the particle shape is oblate spheroidal and when the particles are coated with a dielectric substance such as ice; it will shift to shorter wavelengths when the particle shape is prolate spheroidal (Gilra gilra (1972); Hecht hecht (1981); Draine & Malhotra draine+malhotra (1993)). The observed lack of correlation between the position and FWHM of the peak therefore presents a challenge to the hypothesis that graphite particles originate this peak. In an extensive investigation, Draine & Malhotra (draine+malhotra (1993)) conclude that if graphite particles are the carriers of the 2175 Å peak, a variation in their optical properties ought to be present. This may be the result of varying amounts of impurities, variations in crystallinity, or changes in its electronic structure due to surface effects. Clustering effects have also been considered; Rouleau et al. (rouleau (1997)) have shown that compact clusters of graphitic spheres satisfy the criterion that the position of the peak remain stable while the width varies.

To investigate the clustering effect further we here compute and analyse the extinction of different polycrystalline graphitic clusters. We have chosen clusters ranging from small to large, and that are either sparse or compact. In this way we intend to evaluate how the extinction is influenced by the structure. We focus on clusters consisting of 4, 7, 8, 27, 32, 49, 108 and 343 touching polycrystalline graphitic spheres. The extinction of the clusters is calculated by the use of a rigorous (GA) method (Gérardy & Ausloos GA (1982)) as well as by two DDA methods – MarCODES (Markel markel (1998)) and DDSCAT (Draine & Flatau draine+flatau (1994)) – to test how well these approximations perform when applied to clusters with different morphology. Recently Xu & Gustafson (xu+gustafson (1999)) compared light-scattering calculations by a rigorous solution similar to the GA method, with DDSCAT for two identical polystyrene spheres in contact. They found that DDSCAT worked reasonably well on small volume structures while its validity is challenged on large structures. We find a tendency for both the DDA codes to overestimate significantly the extinction for small (N $`<10`$) fractal clusters and to a much larger extent than for the comparable compact clusters.

By studying clusters consisting of polycrystalline graphitic particles we deal with the anisotropy of graphite in a different way than usual (e.g. Draine draine88 (1988)) and it turns out that this has an influence on the comparison of the calculated extinction with the observed interstellar extinction.

This paper is structured as follows. Section 2 deals with the anisotropy of graphite. Section 3 introduces the fractal and compact clusters studied in this work. Section 4 describes the computational methods under scrutiny; the exact GA solution (Gérardy & Ausloos GA (1982)) and the two different DDA methods, developed by Draine & Flatau (DDSCAT, draine+flatau (1994)) and Markel (MarCoDES, markel (1998)). Section 5 presents and dicusses our results for clusters from one up to 343 particles, in the wavelength range $`0.03100\mu `$m. Section 6 discusses the implications of our results in the interpretation of the interstellar extinction curve. Section 7 presents the conclusions.

## 2 Graphite

Theoretical computation of absorption and scattering by graphitic particles is difficult because graphite is a semi-metal with high anisotropy. Graphite can be characterised by two different dielectric functions, $`ϵ_{}`$ and $`ϵ_{}`$, corresponding to the electric-field vector $`𝐄`$ being perpendicular ($`ϵ_{}`$) and parallel ($`ϵ_{}`$) to the symmetry axis of the crystal ($`c`$-axis), which is perpendicular to the basal plane. It is a lot easier to experimentally determine $`ϵ_{}`$ than $`ϵ_{}`$, because graphite cleaves readily along the basal plane and hence reflectivity measurements can be made with normally incident light, whereas it is very difficult to prepare suitable optical surfaces parallel to the $`c`$-axis.

In this work we deal with the anisotropy of graphite by assuming that in all our clusters, each individual particle is polycrystalline having a dielectric function $`ϵ_{\mathrm{ave}}`$ given by the arithmetic average of $`ϵ_{}`$ and $`ϵ_{}`$, namely $`ϵ_{\mathrm{ave}}=\frac{1}{3}ϵ_{}+\frac{2}{3}ϵ_{}`$. For a polycrystal, this arithmetic average is in fact an attainable upper bound for its dielectric function (Avellaneda et al. avellaneda (1988)). We use the dielectric functions $`ϵ_{}`$ and $`ϵ_{}`$ derived by Draine & Lee (draine+lee (1984)) covering the region from the far-IR to the far-UV. The data agree well in the 0.03$``$62 $`\mu `$m region with those published by Borghesi & Guizzetti (borghesi+guizzetti (1991)).

In constrast, the usual “1/3$``$2/3” approximation, treats individual particles as mono-crystalline - 1/3 of the cluster particles are assumed to have dielectric function $`ϵ_{}`$ and the remaining 2/3 to have dielectric function $`ϵ_{}`$. This approximation has been shown by Draine (draine88 (1988)) and Draine & Malhotra (draine+malhotra (1993)) to have a surprisingly good accuracy for graphite grains with radii $`200`$ Å. However, considering grain formation and grain growth in stellar environments (Sedlmayr sedlmayr (1994)), mono-crystalline particles do not seem to be as valid an assumption as polycrystalline ones. Also, investigations of presolar nano-diamond grains from meteorites - direct specimens of surviving physical material formed in past stellar environments which can be quantitatively analysed in the laboratory - show a clear tendency towards being polycrystalline (Daulton et al. daulton (1994); Phelps phelps (1999)).

A particle in vacuum will show a Lorenz-Mie absorption peak whenever the real part of its dielectric function, $`\mathrm{}(ϵ)`$, satisfies $`\mathrm{}(ϵ)=2`$. Looking at the average dielectric function of graphite shown in Fig. 1, absorption peaks should occur at about $`0.220\mu `$m and $`14.6\mu `$m. The latter should be much more damped due to the higher value of the imaginary part of $`ϵ`$.

## 3 Structure of the clusters

We consider three-dimensional clusters of identical touching spherical particles arranged in three different geometries: fractal, simple cubic, and face-centred cubic. The structures do not have shapes expected to be found in space, but will provide us with boundary conditions for the problem of calculating the extinction from clusters of grains of different morphology. Table 1 lists all the clusters used in this work.

The fractal clusters, frac7, frac49, and frac343, are obtained from the first three stages of the recursive construction of the snowflake fractal. This construction can be summarised as follows: a seed particle is put at the origin of the coordinate system. In the first step a generator is built by symmetrically gluing to the seed particle six copies of itself along the x, y, and z axes. Next, each particle in the first configuration, the generator, is substituted by the whole generator itself. In the next steps the same rule is applied: each particle is replaced by the generator. Fig. 2 shows this procedure up to the second stage. For a fractal cluster its dimension $`D`$ is defined by

$$M(r)=M_0\left[\frac{r}{r_0}\right]^D,$$
(1)
where $`M(r)`$ is the mass of material contained within a sphere of radius $`r`$. For a solid grain of constant density $`D=3`$, whereas for a fractal grain $`D<3`$ (Mandelbrot mandelbrot (1983)); in particular, for the snowflake fractal $`D=\mathrm{ln7}/\mathrm{ln3}=1.77`$ (Vicsek vicsek (1983)).

Although such a deterministic structure (Fig. 2) is not expected to occur in nature, its fractal dimension is close to that of more realistic random cluster-cluster aggregation models. In particular, small particles in space may move in straight, ballistic trajectories and form larger aggregates upon collisions. Numerical simulation of this process yields a value of around 1.9 for the fractal dimension of the resulting aggregates (Meakin meakin (1988); Botet and Jullien botet (1988); Meakin & Jullien meakin+jullien (1988)); this value has also been confirmed by experimental results (Wurm & Blum wurm (1998)). Moreover, optical properties of fractal clusters with $`D<2`$ are predicted to be significantly different from those with $`D>2`$ (Berry & Percival berry (1986)); hence we consider important to use a fractal with a realistic dimension in our computations. The open fractal structure presents a challenge for modelling due to the high porosity and large surface area.

As a contrast to the fractal structure, compact crystalline structures are studied, namely, face-centred cubic and simple cubic, see e.g. Kittel (kittel (1986)) for a discussion on crystal structures. All the clusters listed in Table 1 are symmetric except the fcc49111For such an fcc cluster 48 particles would form a $`3\times 2\times 2`$ unit-cell, while 49 particles can only be fitted into a $`3\times 2\times 3`$ unit-cell.. The non-fractal clusters are a challenge to model since they are far from being spherical.

## 4 The computational methods

### 4.1 Gérardy and Ausloos theory

A rigorous and complete solution to the multi-sphere light scattering problem has been given by Gérardy & Ausloos (GA) (GA4 (1980); GA (1982); GA2 (1983); GA3 (1984)) as an extension of the Mie-Ruppin theory (Mie mie (1908); Ruppin ruppin (1975)). It is based on the exact solution of Maxwell’s equations for arbitrary cluster geometries, polarisation and incidence direction of the light. This is done by expanding the various fields involved in terms of the vector spherical harmonics (VSH). The usual boundary conditions are extended to take into account the possible existence of longitudinal plasmons222A collective excitation for quantized plasma oscillations, in which the free electrons in a metal are treated as a plasma. in the spheres. High-order multipolar electric and magnetic interaction effects are included. We summarise their work as used throughout the paper.

We consider a cluster of $`N`$ homogeneous spheres of radius $`R`$ and dielectric function $`ϵ(\omega )`$, embedded in a matrix of dielectric constant $`ϵ_M`$ and submitted to a plane polarised time harmonic electromagnetic field. The total scattered field from the cluster is represented as a superposition of individual fields scattered from each sphere. The electromagnetic field impinging on each sphere consists of the external incident wave and the waves scattered by the other spheres. For any sphere, the incident, internal and scattered fields are expressed in VSH centred at the sphere origin. The boundary conditions on its surface are solved by transforming all relevant field expansions into the sphere coordinate system, yielding the following system of linear coupled equations for the expansion coefficients $`c_\mu `$ and $`d_\mu `$ of the scattered field (Gérardy & Ausloos GA (1982)):

$`c_\mu `$ $`=\mathrm{\Gamma }_q\{a_{0}^{}{}_{\mu }{}^{}+{\displaystyle \underset{\nu }{\overset{^{}}{}}}(J_{\mu \nu }c_\nu +C_{\mu \nu }d_\nu )\},`$ (2)
$`d_\mu `$ $`=\mathrm{\Delta }_q\{b_{0}^{}{}_{\mu }{}^{}+{\displaystyle \underset{\nu }{\overset{^{}}{}}}(C_{\mu \nu }c_\nu +J_{\mu \nu }d_\nu )\},`$ (3)
where $`\mu =(q,p,i)`$ and $`\nu =(n,m,j)`$. The indices $`q`$ and $`n`$ denote the polar order, $`p=q,\mathrm{},q,m=n,\mathrm{},n`$, and the indices $`i,j`$ number the particles in the cluster. Moreover, $`a_{0}^{}{}_{\mu }{}^{}`$ and $`b_{0}^{}{}_{\mu }{}^{}`$ are the coefficients of the expansion of the external incident wave in VSH in the coordinate frame of the $`i`$ th sphere (Gérardy & Ausloos GA (1982)). The multi-polar electric and magnetic susceptibilities of a sphere (Gérardy & Ausloos GA (1982)) are denoted by $`\mathrm{\Delta }_q`$ and $`\mathrm{\Gamma }_q`$ respectively. The coefficients $`J_{\mu \nu }`$ and $`C_{\mu \nu }`$ describe the transformation of the VSH from a frame centred on particle $`j`$ to another centred on particle $`i`$. Analytical expressions for these coefficients have been given by Gérardy and Ausloos (GA (1982)). The primes in the sums in Eqs. (2) and (3) indicate that terms with $`j=i`$ are omitted. The solution to this system of equations can be written in matrix form as

$$\left[\begin{array}{c}𝐜\\ 𝐝\end{array}\right]=\left[\begin{array}{cc}𝐓_{N11}& 𝐓_{N12}\\ 𝐓_{N21}& 𝐓_{N22}\end{array}\right]\left[\begin{array}{c}𝐚_\mathrm{𝟎}\\ 𝐛_\mathrm{𝟎}\end{array}\right],$$
(4)
where $`𝐓_N`$ is the T-matrix (Waterman waterman (1971); Tsang et al. tsang (1985); Mishchenko et al. 2000b ) of the cluster with components

$`𝐓_{N11}`$ $`=`$ $`((𝚪^1𝐉)𝐂(𝚫^1𝐉)^1𝐂)^1`$ (5)
$`𝐓_{N12}`$ $`=`$ $`𝐓_{N11}𝐂(𝚫^1𝐉)^1`$
$`𝐓_{N21}`$ $`=`$ $`(𝚫^1𝐉)^1𝐂𝐓_{N11}`$
$`𝐓_{N22}`$ $`=`$ $`(𝚫^1𝐉)^1(𝐂𝐓_{N12}+𝐈),`$
and $`𝐈`$ is the identity matrix. When $`N=1`$, Eq. (4) is just the Mie-Ruppin result (Mie mie (1908); Ruppin ruppin (1975))

$$\left[\begin{array}{c}𝐜\\ 𝐝\end{array}\right]=\left[\begin{array}{cc}𝚪& \mathrm{𝟎}\\ \mathrm{𝟎}& 𝚫\end{array}\right]\left[\begin{array}{c}𝐚_\mathrm{𝟎}\\ 𝐛_\mathrm{𝟎}\end{array}\right].$$
(6)
The extinction efficiency, $`Q_{ext}`$, i.e. the extinction cross-section in units of the total geometrical cross-section can then be obtained from Gérardy & Ausloos (GA (1982))

$$Q_{ext}=\frac{1}{\pi Nk^2R^2}\mathrm{}(𝐚_0^T𝐜+𝐛_0^T𝐝),$$
(7)
where $`k=2\pi /\lambda `$ and $`\lambda `$ is the wavelength in the matrix. The extinction cross section per unit volume is then computed as

$$\frac{C}{V}=\frac{3}{4R}Q_{ext}.$$
(8)
Limiting $`q`$ and $`n`$ in Eqs. (2) and (3) to integers less than or equal to $`L`$, we obtain a system of $`2NL(L+2)`$ equations whose solution, Eq. (4), is the $`2^L`$-polar approximation to the electromagnetic response of the cluster; $`𝐚_0,𝐛_0,𝐜\text{ and }𝐝`$ are complex vectors of dimension $`NL(L+2)`$, while the $`𝐓_{Nij},i,j=1,2`$, are complex square matrices of dimension $`NL(L+2)`$. In this case Eqs. (7) and (8) determine the $`2^L`$-polar approximation to the cluster extinction efficiency and the cluster extinction per unit volume respectively.

The extension of this theory to treat the case of absorbing embedding media is given in (Lebedev & Stenzel lebedev1 (1999); Lebedev et al. lebedev2 (1999)).

### 4.2 The DDA method

The discrete dipole approximation (DDA) - also known as the coupled dipole approximation - method is one of several discretisation methods (e.g. Draine draine88 (1988); Hage & Greenberg hage+greenberg (1990)) for solving scattering problems in the presence of a target with arbitrary geometry. The discretisation of the integral form of Maxwell’s equations is usually done by the method of moments (Harrington harrington (1968)). Purcell & Pennypacker (purcell (1973)) were the first to apply this method to astrophysical problems; since then, the DDA method has been improved greatly by Draine (draine88 (1988)), Goodman et al. (goodman (1991)), Draine & Goodman (draine+goodman (1993)), Draine & Flatau (draine+flatau (1994)), Markel (markel (1998)), and Draine (draine00 (2000)). The DDA method has gained popularity among scientists due to its clarity in physical principle and the FORTRAN implementation which have been made publicly available by Draine & Flatau (DDSCAT package, Draine & Flatau draine+flatau (1994)) and by Markel (MarCoDES, Markel markel (1998)).

When considering the problem of scattering and absorption of linearly polarised light

$$E_0=e_0\mathrm{exp}\left(ikr\right)$$
(9)
by an isotropic graphitic grain, then within the concept of the DDA method, the grain is replaced by a set of discrete elements of volume $`V_i`$ with relative dielectric constant $`ϵ_i`$ and dipole moments $`d_i=d\left(r_i\right)`$, $`i=1,\mathrm{},N,`$ whose coordinates are specified by vectors $`r_i`$. The equations for the dipole moments can be written using simple considerations based on the concept of the exciting field, which is equal to the sum of the incident wave and the fields of the rest of the dipoles in a given point

$$𝐝_i\left(r_i\right)=\alpha _i\left[E_0\left(r_i\right)+k^3\underset{ji}{}𝐆_{ij}𝐝_j\left(r_j\right)\right]$$
(10)
where the dipole scattering tensor $`𝐆_{ij}`$ has the following form in dyad notations:

$`𝐆_{ij}`$ $`=`$ $`\left[A\left(kr\right)\delta _{ij}+B\left(kr\right){\displaystyle \frac{r_ir_j}{r^2}}\right]`$ (11)
$`A\left(x\right)`$ $`=`$ $`\left[x^1+ix^2x^3\right]\mathrm{exp}\left(ix\right)`$
$`B\left(x\right)`$ $`=`$ $`\left[x^13ix^2+3x^3\right]\mathrm{exp}\left(ix\right)`$
$`𝐫`$ $``$ $`𝐫_{ij}=𝐫_i𝐫_j.`$
The solution of the set of Eqs. (10) and (11) yields all the basic optical characteristics of a particle such as the integrated extinction $`\left(Q_{\mathrm{ext}}\right)`$ and absorption $`\left(Q_{\mathrm{abs}}\right)`$:

$`Q_{\mathrm{ext}}`$ $`=`$ $`4\pi k\mathrm{}\left[{\displaystyle \underset{i}{}}\left(e_0d_i\right)\mathrm{exp}(ik_0r_i)\right]`$ (12)
$`Q_{\mathrm{abs}}`$ $`=`$ $`4\pi k{\displaystyle \underset{i}{}}\eta _i\left|d_i\right|^2,\eta _i={\displaystyle \frac{4\pi \mathrm{}\left(ϵ_i\right)}{V_i\left|ϵ_{i1}\right|^2}},`$
and scattering $`\left(Q_{\mathrm{sca}}=Q_{\mathrm{ext}}Q_{\mathrm{abs}}\right)`$ efficiencies. Here $`\mathrm{}`$ means the imaginary part of the expression in the argument.

Once the location and polarisability of the points are specified, calculations of the scattering and absorption of light by the array of polarisable points can be carried out to in principle whatever accuracy is required. The limiting factor is the capacity of computing resources.

#### 4.2.1 DDSCAT

In this work we use the DDSCAT code version 5a10 (Draine & Flatau draine+flatau (1994); Draine & Flatau draine+flatau00 (2000)). This version contains a new shape option where a target can be defined as the union of the volumes of an arbitrary number of spheres. In DDSCAT the considered grain/cluster is replaced by a cubic array of point dipoles. The cubic array has numerical advantages because the conjugate gradient method can be efficiently applied to solve the matrix equation describing the dipole interactions (Goodman et al. goodman (1991)). When knowing the dipole strength of each cell in the particle it is possible to compute the optical properties of arbitrary dust configurations.

There are three criteria for validity of DDSCAT:

1. The wave phase shift $`\rho =|m|kd`$ ($`m=\sqrt{ϵ}`$ being the complex refractive index of the target material) over the distance $`d`$ between neighbouring dipoles should be less than 1 for calculations of total cross sections and less than 0.5 for phase function calculations.
2. $`d`$ must be small enough to describe the object shape satisfactorily.
3. The refractive index $`m`$ must fulfill $`|m|<2`$.

A comparison study by Draine & Flatau (draine+flatau (1994)) shows that scattering and absorption cross sections can be calculated with DDSCAT to accuracies of a few percent, provided that the criteria elaborated above are satisfied.

Xing & Hanner (xing+hanner (1997)) state that the $`Q`$ efficiency factor is generally much less sensitive to those criteria, and that it is usually sufficient to let the value of $`\rho `$ be around 1 or possible bigger. Xu & Gustafson (xu+gustafson (1999)) show through a comparison between the exact solution for two spheres in contact and DDSCAT that the criterion set up by Draine & Flatau (draine+flatau (1994)) is not sufficient to ensure high accuracy if the particles are large and strongly interacting (high refractive index). They, therefore, recommend to use a validity criteria of $`\rho <0.3`$ when calculating total cross sections.

The typical number of dipoles needed to obtain a reliable computational result using the DDSCAT code can be determined by calculating the minimum number of dipoles needed per particle. When a particle is represented by a 3-dimensional array of $`N`$ dipoles, its volume is $`Nd^3`$, which must be equal to $`4\pi R^3/3`$,

$$N=\frac{4\pi }{3}\left(\frac{R}{d}\right)^3=\frac{4\pi }{3}\left(\frac{2\pi R|m|}{\rho \lambda }\right)^31039\left(\frac{R|m|}{\rho \lambda }\right)^3$$
(13)
since $`d`$ is related to the wave phase shift $`\rho `$ by $`d=\rho /(|m|k)`$ (Draine & Flatau draine+flatau00 (2000); Xing & Hanner xing+hanner (1997)). The number of dipoles needed to satisfactorily describe the shape of the object considered might be larger than indicated by Eq. (13) if the shape is “challenging”. It is important to remember that Eq. (13) only fulfills the first of the three criteria set by Draine & Flatau (draine+flatau00 (2000)).

For materials with large refractive index ($`|m|>2`$), Draine & Goodman (draine+goodman (1993)) show that especially the absorption is overestimated by DDA. The limitation in DDSCAT is set by the use of the lattice dispersion relation (LDR) for electromagnetic waves propagating on an infinite cubic lattice of point dipoles of polarisability $`\alpha _i`$ and spacing $`d`$. According to Draine & Goodman (draine+goodman (1993)), the LDR prescription for $`\alpha _i`$ gives fair accuracy for scattering but poorer results for absorption. When $`|m1|`$ is large the continuum material is effective at screening the external field: in the limit $`|m1|\mathrm{}`$ the internal field generated by the polarisation would exactly cancel the incident field, so that the continuum material in the interior of the target would be subjected to zero field. In the case of a discrete dipole array, the dipoles in the interior will also be effectively shielded, while the dipoles located on the target surface are not fully shielded and, as a result, absorb energy from the external field at an excessive rate. This error can be reduced to any desired level by increasing the number $`N`$ of dipoles, thereby minimising the fraction $`N^{1/3}`$ of the dipoles which are at surface sites, but very large values of $`N`$ are required when $`|m|`$ is large (Draine & Goodman draine+goodman (1993)).

Graphite is characterised by having a high refractive index. This means that the criterion $`|m|<2`$ is only fulfilled shortward of 0.072 $`\mu `$m ($`m=n+ik=0.58+i1.42`$) and relaxing the criterion a little (to account for the variability of $`m`$ that occurs) gives an upper limit of 0.216 $`\mu `$m ($`m=0.66+i1.35`$).

#### 4.2.2 MarCoDES

Another efficient code based on DDA is the Markel Coupled Dipole Equation Solver (MarCoDES, Markel markel (1998)). This code is designed to solve the coupled Eq. (10) for an arbitrary cluster of $`point`$ dipoles using either conjugate gradient method (iterative method) or the LU expansion (direct method). The program is in principle applicable to clusters of small particles with arbitrary geometry, but is most computationally efficient for sparse clusters (i.e. when the volume fraction is very low) with significant number of particles $`(10^310^4)`$. Contrary to DDSCAT the program does not use the Fast Fourier Transformation (FFT) because this might significantly decrease the computational performance for clusters with a low volume filling factor. When the volume filling fraction is close to unity, algorithms utilising FFT will be much faster. In the framework of the MarCoDES a spherical grain is considered to be equivalent to the single elementary dipole, this corresponds to $`N=1`$ in the DDSCAT case. The dimensionality of the coordinates of particles in MarCoDES require a special consideration. By replacing real particles by point dipoles located at their centres we significantly underestimate the strength of their interaction. In order to correct the interaction strength, the author of MarCoDES introduces geometrical intersection of particles. All coordinates are defined in terms of the distance between neighbouring dipoles $`d`$, which is given by $`d=\left(4\pi /3\right)^{1/3}R`$. So, for example, if two particles have radii 10 nm, then the distance between the dipoles is 16.12 nm. This suggested phenomenological procedure allows MarCoDES to be more accurate than the usual single dipole approximation since the intersection produces some analogy of including higher multipole interactions between particles. Application of MarCoDES to fractal sparse clusters allows, in principle, to introduce renormalisation on the sphere radii and number of particles depending on the fractal dimension of the cluster (Markel et al. markel00 (2000)), which will produce even further corrections to the strength. In this paper we use the simple particle intersection model given in the MarCoDES. The fact that the program only uses a single dipole for each particle in the cluster have significant benefits in computation efficiency when compared to other multi-polar approaches such as the GA method and DDSCAT.

## 5 Results

We present the results of the computation of the extinction of infrared, visible, and ultraviolet light by clusters of graphitic spheres of radii 10 nm and 50 nm, containing as few as 1 sphere to as many as 343 spheres. Specifically, our results concern the wavelength region from 0.1 to 100 $`\mu `$m. For convenience, we shall, from now onwards, call a cluster small if it has less than 10 particles, medium-size if it has between 25 and 50 particles, and large if it has more than 100 particles. With this convention, we shall first discuss the convergence of the GA computations, followed by a comparison of the three computational methods (GA, MarCoDES, and DDSCAT) when applied to the calculation of the extinction of a single polycrystalline graphitic sphere. Then we shall present in succession our results for small, medium-size, and large clusters, and end the section with an overall assessment of the computational methods in light of the results presented.

### 5.1 Convergence of GA method

As we stated in Sect. 4.1, with the GA method we can compute the extinction of a cluster in the $`2^L`$-polar approximation. We find that the smallest L needed for the convergence of the extinction of our clusters will be different for different regions of the optical spectrum. For small clusters, our full converged results show that in the UV-visible it is sufficient to use L = 7 for open clusters, and L = 9 for compact ones, to compute the extinction with an accuracy of 1%; whereas for longer wavelengths, at least L = 11 is required to ensure the same accuracy; this is so because in this region graphite has a metallic-like behaviour (high $`|m|`$), and, thus, higher order multipoles have to be included in the computation of the extinction to get the same accuracy as in the UV-visible range. We note that in the UV-visible range, by accepting an accuracy of 5% in the computation of the extinction, we can use L = 5 for open clusters and L = 7 for compact ones; we expect this to hold for clusters of up to a few tens of particles.

In column 2 of Table 2, we give the cut-off polar order L used in the calculation of the extinction of all our clusters. Full convergence was achieved only for small clusters; for the others, L indicates the maximum polar order obtained with our computer resources. In general, L = 11 was used for small clusters, L from 6 to 7 for medium-size clusters, and L from 2 to 3 for large clusters.

We show in Fig. 3 the convergence of the extinction of the frac7 cluster in the wavelength range $`0.03100\mu `$m. The inset shows the 2175 Å peak; around the peak maximum L = 3 already gives results accurate within 1%.

### 5.2 Single sphere

To set up a comparison baseline, we compute the extinction of af single sphere of radius 10 nm using the two DDA codes (DDSCAT and MarCoDES) and the GA solution. Since the GA theory extends that of Mie (Sect. 4.1), the GA and Mie solution coincide for a single sphere; this solution will be taken as a reference for the comparison.

It can be seen from Fig. 4 that DDSCAT is almost indistinguishable from the Mie solution for $`\lambda 1.0\mu `$m while MarCoDES is so for $`\lambda 0.6\mu `$m. The less good fit by DDSCAT for $`\lambda >1.0\mu `$m is most likely a consequence of the refractive index of graphite having a modulus larger than 2 for these wavelengths, and of the 24 464 dipoles being inadequate to account for this.

### 5.3 Small ($`N<10`$) clusters

As before, we take the fully converged GA solutions as a reference for comparison. Fig. 5 shows the extinction of the frac7 cluster as computed with the three methods, while Fig. 6 shows similar information for the sc8 cluster. For the frac7 cluster (Fig. 5) MarCoDES underestimates the extinction for $`0.2\mu `$m $`<\lambda <0.25\mu `$m and systematically overestimates it for $`\lambda >0.26\mu `$m. DDSCAT333Calculated with a $`36\times 36\times 36`$ dipole grid, which provides 6838 dipoles ($``$ 977 dipoles per particle). slightly underestimates the extinction for $`\lambda <0.23\mu `$m and overestimates it for $`\lambda >0.23\mu `$m. In the range $`0.26\mu `$m $`<\lambda <4\mu `$m the excess extinction from DDSCAT is smaller than for the MarCoDES calculations while for $`\lambda >4\mu `$m DDSCAT gives a systematic relative increase in the excess extinction with wavelength which is much larger than the almost constant factor, $`1.5`$, seen for MarCoDES.

For the sc8 cluster (Fig. 6) the trend is different. DDSCAT444Calculated with a $`32\times 32\times 32`$ dipole grid, which provides 16 824 dipoles ($``$ 2103 dipoles per particle). overestimates the extinction for $`\lambda <0.3\mu `$m, underestimates it for $`0.3\mu `$m $`<\lambda <11\mu `$m and overestimates it again for $`\lambda >11\mu `$m. MarCoDES also overestimates the extinction for $`\lambda <0.26\mu `$m but then systematically underestimates it to a much larger extent than DDSCAT for $`\lambda >0.26\mu `$m.

We have also investigated the effect of particle-size on a cluster’s extinction using the fully converged results of the GA solution. For this, in addition to the extinction of the frac7 cluster of spheres of radii 10 nm, we have also computed the extinction of a frac7 cluster of spheres of radii 50 nm; which was fully converged at L = 11 with a precision of 1.2%. Fig. 7 shows the results; what strikes the eye immediately is the high increase in the extinction over almost the entire wavelength range investigated, upon increasing the radii of the spheres. Still, for wavelengths $`\lambda >1.1\mu `$m the shape of the extinction of both clusters look roughly the same. More important, however, are the differences observed in the wavelength range $`\lambda <1.1\mu `$, where the extinction of the cluster with the larger spheres exhibits peaks at $`0.11\mu `$m and $`0.4\mu `$m in addition to the 2175 Å peak, which is the only peak displayed by the extinction of the cluster with the smaller spheres. All these observations agree with the finding of Kimura (2001) that the light scattering of fractal clusters depends strongly on the size parameter of the monomer.

### 5.4 Medium-size ($`25<N<50`$) clusters

For the frac49 cluster, shown in Fig. 2, the extinction cross section obtained by the three methods is displayed in Fig. 8. For this structure neither DDSCAT555Calculated with a $`36\times 36\times 36`$ dipole grid, which provides 1722 dipoles ($``$ 35 dipoles per particle). nor MarCoDES gets close to the solution suggested from the GA calculations with L = 6. The GA solution is not fully converged at long wavelengths due to computational limitations; however, in the UV-visible region convergence is ensured at least with an accuracy of 5%. Both MarCoDES and DDSCAT greatly overestimates the extinction cross section for this cluster at longer wavelengths.

As a comparison to the frac49 cluster we have carried out calculations for two symmetric clusters, sc27 and fcc32, and an asymmetric one, fcc49.

For the sc27 cluster (Fig. 9) the GA solution was truncated at L = 7. MarCoDES underestimates the extinction for $`0.17\mu `$m $`<\lambda <0.22\mu `$m and $`\lambda >0.29\mu `$m. DDSCAT666Calculated with a $`32\times 32\times 32`$ dipole grid, this provides 16 783 dipoles ($``$ 622 dipoles per particle). also slightly underestimates the extinction for $`0.17\mu `$m $`<\lambda <0.22\mu `$m and $`0.26\mu `$m $`<\lambda <0.4\mu `$m, but systematically overestimates it for $`\lambda >1\mu `$m.

We found that the version of MarCoDES tested here can not calculate a face-center cubic structure of touching particles because in this case the lattice cells representing neighbouring particles will touch only at the corners, giving as a consequence a spectrum corresponding to non-touching particles. The extinction of fcc32 and fcc49 clusters were therefore only calculated with the GA theory and DDSCAT. The results for fcc32 are displayed in Fig. 10; the results for the fcc49 cluster look similar. The GA calculations were truncated at L = 7 for fcc32 and at L = 6 for fcc49. For the DDSCAT calculations a $`32\times 32\times 32`$ dipole grid was used; for the fcc32 cluster this provides 17 942 dipoles ($``$ 561 dipoles per particle). For fcc49 a $`48\times 32\times 40`$ dipole grid was used; this provides 29 849 dipoles ($``$ 609 dipoles per particle). As can be seen in Fig 10, DDSCAT systematically overestimates the extinction especially at long wavelengths. The overestimation of the extinction cross section as calculated by DDSCAT is a general tendency for all the clusters we have studied. Also DDSCAT does not reproduce the absorption hump expected for graphite around $`1015\mu `$m, which is reproduced in the GA and MarCoDES calculations and also in Lorenz-Mie calculations by Draine & Li (draine+li (2001)).

### 5.5 Large ($`N>100`$) clusters

The GA solution could only be calculated to L = 3 for the fcc108 cluster and to L = 2 for the frac343 cluster due to computational limitations. The extinction of the fcc108 and frac343 clusters looks very similar to that of the fcc32 and frac49 clusters shown in Figs. 10 and 8, respectively. For the frac343 cluster the whole spectrum is shifted to longer wavelengths by about 0.002 $`\mu `$m and for the fcc108 cluster the whole spectrum is shifted 0.01 $`\mu `$m to longer wavelengths compared to the frac49 and fcc32 cluster, respectively. It was not possible to calculate these large clusters with DDSCAT.

### 5.6 Discussion

It is seen that both DDA codes overestimate the extinction of the fractal structures at long wavelengths. A likely explanation is that the DDA model considers only point dipole interactions, whereas the GA solution matches the boundary conditions on the surface of the spherical particle and solves the problem of scattered and internal fields at the same time. Also, the large surface area of fractal clusters contributes to the complexity of the scattered field, which the GA method handles well because it takes into account the surface topology through the boundary conditions; the DDA codes, on the other hand, perform best when the surface to volume ratio is low. This interpretation is also supported by the calculations for the small compact sc8 cluster, see Fig. 6, where the DDA approximations perform reasonably well. For larger compact clusters, however, DDSCAT does not give good results at long wavelengths. This is due to computational limitations; the number of dipoles per particle that we could handle becomes lower as the cluster size increases. It is interesting to note that MarCoDES gives better results for compact than for fractal cluters of intermediate size. According to Markel et al. (markel00 (2000)), the performance of MarCoDES can be improved by altering the intersection parameter which determines if the particles touch or overlap, but we have not investigated that here. Nevertheless, the fact that MarCoDES in some cases overestimates the extinction, as compared to GA theory, and in other cases underestimates it, suggests that the determination of the rather arbitrary optimal intersection parameter is a very complex problem indeed.

We next compare the extinction of fractal and compact clusters of intermediate size. In Fig. 11, the GA calculation for frac49 is compared to that for fcc49 and fcc32. It is seen that the extinction of fractal and compact clusters are of the same order of magnitude. This is in contradiction to the DDA calculations which give a much higher extinction for fractal clusters (Figs. 5 \- 11). We, therefore, recommend that great care should be taken when using these DDA codes at longer wavelengths for fractal clusters with low fractal dimension consisting of materials with high refractive index ($`m`$).

Concerning the effect of particle size on the extinction cross section, in classical electromagnetic theory the extinction is independent of cluster size in the long wavelength limit, i.e. when the size is much smaller than the wavelength. As seen in Fig. 7, clusters consisting of particles of radii 50 nm are clearly outside of this limit. Decreasing the particle radius below 10 nm, however, will give, in a classical theory, at most minor effects on the extinction of single particles and small clusters. Quantum mechanical effects, on the other hand, may influence the width and the position of the absorption peak for such small particles. A detailed study of these particle size effects on the extinction cross section is outside the scope of the present work.

## 6 The 2175 Å extinction peak

The main observational constraints concerning the 2175 Å peak are according to Fitzpatrick & Massa (fitzpatrick+massa86 (1986)) and Mathis (mathis (1994)):

1. The peak position is remarkably constant, $`\nu _{\mathrm{max}}=4.6\pm 0.04\mu `$m-1, where $`\nu =1/\lambda `$.
2. The peak width has a much larger range of variations of $`1.0\pm 0.25\mu `$m-1 (i.e. $`25\%`$).
3. The variation in peak position and the width ($`\gamma `$) are uncorrelated, except for the widest bumps ($`\gamma 1.2\mu `$m-1) where a systematic shift to larger wavenumber is observed (Cardelli & Savage cardelli+savage (1988); Cardelli & Clayton cardelli+clayton (1991)). These later observations all have lines of sight passing through dark, dense regions.

Table 2 lists the peak position and width for the different cluster structures, as computed with the GA method. The polar order at which the GA calculation was truncated is also indicated. The width was determined as the FWHM. For the asymmetric compact cluster, we have investigated the importance of the orientation of the cluster and it was found that the peak position was shifted less than 0.1 $`\mu `$m-1 while the width of the feature was not affected.

It was shown in a study by Rouleau et al. (rouleau (1997)), that small compact clusters qualitatively satisfy the observational constraints on the 2175 Å peak, except that the peak position falls at the wrong wavelength. This result is not confirmed from our findings, since the compact clusters we have studied do not have a stable peak position. The clusters considered by Rouleau et al. (rouleau (1997)) all had peak positions at higher wavenumber than the 2175 Å peak while the clusters we have considered all have peak position at lower wavenumbers. Rouleau et al. (rouleau (1997)) also used optical constants from Draine & Lee (draine+lee (1984)), but they dealt with the anisotropy of the graphitic particles with the usual “$`1/32/3`$” approximation, see Sect. 2, which assumes the particles to be mono-crystalline. As argued in Sect. 2 assuming the particles to be polycrystalline seems much more reasonable. This implies that variations in crystallinity is indeed an important factor when discussing the 2175 Å peak as also suggested by Draine & Malhotra (draine+malhotra (1993)).

The fractal clusters considered here do not enhance the long wavelength wing and do not introduce additional structure as seen by Rouleau et al. (rouleau (1997)) for elongated and ”fluffy” structures. The frac7 and frac49 clusters do come close (peak position is off by $`0.040.08\mu `$m-1) to the observational constraints. This indicates that small ($`N5100`$) fractal clusters ought to be investigated in more detail to determine if fractal clusters with low fractal dimension have a stable peak position around $`4.6\mu `$m-1 and produce a variable width depending on the number of particles in the cluster.

## 7 Conclusions

We computed the extinction of clusters of polycrystalline graphitic particles in the wavelength range $`0.1`$ to $`100\mu `$m. Computations by the rigorous multipolar theory of Gérardy & Ausloos (GA (1982); GA) were compared to two formulations of the discrete dipole approximation (DDA).

We have compared the extinction of open fractal clusters and compact clusters. The fractal and compact clusters display an extinction at long wavelengths, of the same order of magnitude as computed with the GA method. It seems that the DDA approximations grossly overestimate the long-wavelength extinction of small fractal structures. If this also holds true for larger fractal clusters is not possible for us to say due to computational limitations. The DDA codes should, however, be used with caution for this type of problem.

The GA computations were compared to the observed interstellar extinction curve. We found that results for small and medium-size (less than 50 particles) fractal clusters are in fair agreement with observational constraints, while those of compact sc and fcc clusters are not.

Convergence of the GA computations for small (less than 10 particles) clusters was found by including 11 multipoles. Good approximative results were obtained for medium-size (between 25 and 50 particles) clusters. In most cases we found a substantial differences between the GA theory and the DDA approximations of Draine & Flatau (draine+flatau (1994); DDSCAT) and Markel (markel (1998); MarCoDES).

As predicted by Draine & Goodman (draine+goodman (1993)) the absorption is significantly overestimated by DDSCAT at the wavelengths where the modulus of the refractive index, $`|m|`$, is large, see Figs. 8, 9 and 10. For $`|m|<2`$, DDSCAT performs well and better than the MarCoDES. MarCoDES, however, is computationally much faster than DDSCAT and the GA calculations. With DDSCAT the accuracy is directly dependent on the number of dipoles used in the approximation. Hence, by using more dipoles than we have used here, a better accuracy could be obtained, but then the computational cost would increase. Which DDA code is best to use will therefore depend on the type of problem one wants to address and the accuracy needed.

###### Acknowledgements.

We would like to thank B.T. Draine and V.A. Markel for making their DDA codes available as shareware. VP would like to thank V.A. Markel and ACA would like to thank Björn Davidsson for fruitful discussions. ACA acknowledges support from the Carlsberg Foundation and from NorFA. JS and GN acknowledges support from the Swedish Natural Science Research Council (NFR). VP acknowledges support from the Wenner-Gren Foundation.


format: text

text:
# Constraints Imposed by a High Magnetic Field on Models for the EUV Emission in the Coma Cluster

## 1 INTRODUCTION

EUV emission in the Coma cluster in excess of that produced by the thermal X-ray gas was first reported by Lieu et al. (1996). Further analysis of data on this cluster was carried out by Bowyer and Berghöfer (1998) and Bowyer et al. (1999). Initially, this flux was attributed to thermal emission from ”warm” gas at $`10^6`$ K. Maintenance of a warm intracluster gas is extraordinarily difficult, and on these grounds alone it was generally agreed that a thermal source was untenable. Observations relevant to this issue were obtained with the Hopkins Ultraviolet telescope (Dixon et al. 1996), and FUSE (Dixon et al. 2001a,b). No Far UV line emission from gas at $`10^6`$ K was detected. More recently, observations of several clusters with XMM (Tamura et al. 2001; Kaastra et al. 2001; Peterson et al. 2001) have shown no evidence for a $`10^6`$ temperature gas. All additional clusters examined with XMM also show no evidence for a “warm” $`10^6`$ gas (S. Kahn, private communication). The sum of all these findings seems compelling: a thermal mechanism for the EUV excess can be ruled out. (However, for an alternate point of view see Mittaz et al. 1998; Lieu, Bonamente, & Mittaz 1999; Lieu et al. 1999; Lieu, Bonamente, & Mittaz 2000; and Bonamente, Lieu, & Mittaz 2001a,b.)

Since the source mechanism is not thermal, it must be the product of a non-thermal process. Inverse Compton scattering (ICS) of cosmic rays with the 2.7 K cosmic microwave background was suggested early-on as a possible source mechanism (Hwang 1997; En$`\beta `$lin & Biermann 1998). This is the only non-thermal mechanism that has been suggested as a possibility for the source of this flux.

A number of researchers have explored models of ICS in an attempt to explain the EUV excess in the Coma cluster (Hwang 1997, Bowyer & Berghöfer 1998, En$`\beta `$lin & Biermann 1998, Sarazin 1999, Atoyan & Völk 2000, Brunetti et al. 2001). With one exception (Atoyan & Völk 2000), all of these models require a low ($`<1\mu `$G) field and typically these models require a very low magnetic field ($`<0.5\mu `$G) field.

A hard (25-80 keV) X-ray flux in excess of that produced by the $`10^8`$ K thermal X-ray gas has been detected with observations with $`BeppoSAX`$ (Fusco-Femiano et al. 1999). Brunetti et al. (2001) have attempted to explain both the EUV and hard X-ray excesses in Coma by a single complex model. This model also employed a low magnetic field.

New studies of the magnetic field in clusters of galaxies strongly suggest that magnetic fields in clusters of galaxies are $`5\mu `$G. In this paper we investigate constraints imposed by a high magnetic field on models capable of producing the EUV excess in clusters of galaxies via the ICS process. We also briefly consider the effects of a high magnetic field on attempts to produce both the EUV and the hard X-ray emission by the ICS mechanism with the same population of cosmic ray electrons.

## 2 ASSUMPTIONS AND OVERALL APPROACH

We make the following assumptions.

A.) The magnetic field is $`5\mu `$G. The magnetic field in the Coma cluster has been the subject of considerable study and observational results vary considerably. Estimates based in part on equipartition arguments typically result in low field estimates (Rephaeli 1988, Giovannini et al. 1993). However, recent observational results strongly indicate that cluster magnetic fields are quite large. Clarke et al. (2001) studied sixteen clusters; a substantial number of data points were obtained for each cluster. They found the cluster magnetic fields were all in the range $`47\mu `$G. These results seem compelling, but it could be argued that they are not applicable to the Coma cluster since that cluster was not included in their study. However, the Coma cluster field can be determined by a different procedure albeit for only one location in the cluster. Feretti et al. (1995) studied the rotation measure of the cluster-embedded head-tail radio source NGC 4869 (5C4.81). They found the number of magnetic field reversals through the cluster $`>200`$ indicating B $`>4.9\mu `$G. This measurement is limited by the resolution of the VLA configuration employed and, in principle, is a lower limit for the number of magnetic field reversals through the cluster. It is possible that with even higher resolution the number of magnetic field reversals through the cluster could be larger, with a resulting increase in the estimate of B. However, this estimate is within the 4 to 7 $`\mu `$G range found by Clarke et al. (2001) for all clusters in their sample. In view of these observations, we find the assumption of 5$`\mu `$G for the field in the Coma cluster to be quite reasonable.

B.) A scattering process is operative in the cluster which will result in a sufficiently long path length for the electrons that the deposition of the energy of the electrons occurs within the cluster itself. We note there is overwhelming observational support for this assumption since if the cosmic rays did not deposit their energy within the cluster, they would escape the cluster and produce excess EUV emission over a region which is far larger than the cluster itself.

C.) We assume the evolution of the electron population is given by the equation (Sarazin 1999):

$$\frac{N(\gamma )}{t}=\frac{[b(\gamma )N(\gamma )]}{\gamma }+Q(\gamma )$$
(1)
where $`N(\gamma )d\gamma `$ is the total number of electrons in the range $`\gamma `$ to $`\gamma +d\gamma `$, $`b(\gamma )`$ is the energy loss rate of an electron with energy of $`\gamma `$, $`Q(\gamma )d\gamma `$ is the injection rate of electrons in the energy range $`\gamma `$ to $`\gamma +d\gamma `$.

D.) We assume the injected cosmic-ray electrons follow a power law distribution.

E.) We assume the only energy loss mechanisms are inverse Compton emission, synchrotron radiation, Coulomb ionization and bremsstrahlung.

The energy loss rate for relativistic electrons interacting with the cosmic microwave background is given by Longair (1994):

$$b_{IC}(\gamma )=\frac{4}{3}\frac{\sigma _T}{m_ec}\gamma ^2U_{CMB}=1.37\times 10^{20}\gamma ^2\mathrm{s}^1$$
(2)
where $`\sigma _T`$ is the Thomson cross section and $`U_{CMB}`$ is the energy density of the cosmic microwave background. The loss rate by synchrotron emission is given by Longair (1994):

$$b_{syn}(\gamma )=\frac{4}{3}\frac{\sigma _T}{m_ec}\gamma ^2U_B=1.30\times 10^{21}\gamma ^2(\frac{B}{1\mu G})^2\mathrm{s}^1$$
(3)
where $`B`$ is the magnetic field and $`U_B=B^2/8\pi `$ is the energy density of the magnetic field. The loss rate for Coulomb ionization is given by Sarazin (1999):

$$b_{Coul}(\gamma )1.2\times 10^{12}n_e[1.0+\frac{\mathrm{ln}(\gamma /n_e)}{75}]\mathrm{s}^1$$
(4)
where $`n_e`$ is the thermal electron density in the ICM of the Coma cluster. The loss rate for the Bremsstrahlung is approximately (Sarazin 1999):

$$b_{brem}(\gamma )1.51\times 10^{16}n_e\gamma [\mathrm{ln}(\gamma )+0.36]\mathrm{s}^1$$
(5)
A comparison of the energy loss rates for these four mechanisms is given in Fig.1. This shows that ICS is the dominant energy loss rate mechanism for relativistic electrons in a $`0.2\mu `$G cluster field while synchrotron radiation is the main energy loss rate process for relativistic electrons in $`5\mu `$G cluster field. A diagram of the energy loss rate timescale $`\gamma /(b_{IC}+b_{syn}+b_{Coul}+b_{brem})`$ versus electron energies is shown in Fig.2. Note that there is a decrease of 10$`\%`$ in the energy of the peak energy loss, and, more importantly, there is a decrease in the peak energy loss rate of $`40\%`$ for the case of a $`5\mu `$G field as compared to a $`0.2\mu `$G field.

The observed EUV ICS flux is given by Blumenthal & Gould (1970):

$$f_{IC}(ϵ)=\frac{1}{4\pi D^2}\times \frac{r_o^2}{\mathrm{}^3c^2\pi }\times K_e\times (kT)^{\frac{p+5}{2}}\times F(p)\times ϵ^{\frac{p+1}{2}},$$
(6)
where the symbols have their standard meaning. The energy of the photon observed, $`ϵ`$, has been obtained by folding the effective area of the EUVE Deep Survey Telescope (Sirk et al. 1997) with absorption of EUV emission by the Galactic interstellar medium in the direction of Coma following the prescription of Bowyer et al. (1999). The peak of the resultant telescope effective area is at 80 Å.

The synchrotron radio flux is given by Blumenthal and Gould (1970):

$$f_{syn}(\nu )=\frac{1}{4\pi D^2}\times \frac{4\pi e^3}{m_ec^2}\times K_e\times B^{\frac{p+1}{2}}\times (\frac{3e}{4\pi m_ec})^{\frac{p1}{2}}\times a(p)\times \nu ^{\frac{p+1}{2}}.$$
(7)
We have used an integrated diffuse flux density for Coma C of 640 mJy at 1.4 GHz (Deiss et al. 1997) in combination with a spectral index $`\alpha =(p1)/2=1.16`$ (Bowyer & Berghöfer 1998) as our observational base.

We investigate three scenarios. In the first, we consider a large injection event occurring over a relatively short time which produces the cosmic-ray electron population. In the second scenario, we assume a continuous sequence of small injection events. Finally, we consider a combination of these two: a large injection event followed by continuous smaller events. In Table 1 we show the input parameters for these three models.

In this Table, $`N_o`$ is the normalization factor and $`p`$ is the spectral index of the initial injected rate. $`Q_o`$ refers to the normalization factor and $`p^{}`$ is the spectral index of the continuous injection rate. $`Q_o`$ is related to the injected electron number $`Q`$ by the relation $`Q=Q_o\times \gamma ^p^{}`$. $`n_e`$ is the thermal electron density in the ICM. EI denotes the energy input.

## 3 RESULTS

We first consider the case of a single primary injection event. The results with input parameters from Model 1 are shown in Fig.3. This population will produce the EUV emission as reported by Bowyer et al. (1999) shown as a cross. This population greatly overproduces the observed synchrotron radio emission which is shown as a heavy solid line in the figure. With the passage of time, synchrotron losses degrade the higher energy electrons, and the resultant population after $`1.4\times 10^8`$ yr and after $`2.5\times 10^8`$ yr is shown by line A and by line B respectively. The results shown in this figure demonstrate why it is not possible to produce both the EUV and the radio emission by a population of electrons produced in a single initial event.

We next consider the case of a continuous distribution of smaller events. The results with the input parameters from Model 2 are shown in Fig.4. The dashed line shows the total injected electrons. The thin solid line shows the electron population after $`2\times 10^9`$ yr. The heavy solid line shows the distribution of the electrons that are required to produce the observed synchrotron radio emission and the cross shows the distribution of the electrons producing the observed EUV emission. The results shown in this figure show why no uniform continuous distribution of injection events can reproduce the observational data.

We have independently derived a variety of electron populations using a cluster field of $`0.2\mu `$G; these are all capable of reproducing the EUV and radio observations. This is consistent with results found by a number of previous studies that demonstrated that a continuous distribution of injection events can produce populations that will replicate both the EUV and the radio data in the case of low cluster fields.

We next consider a primary event followed by continuous injections. We use input parameters listed in Models 3 through 11. The population obtained using the parameters in Model 3 after an evolution time of 1 Gyr are shown as a dashed line in Fig.5. We also show the effect of changes in the size of the initial event. These results are shown as line A and employ inputs from Model 4. Finally, we show the effect of changes in $`p`$. Using inputs taken from Model 5, we obtain the population shown as line B.

In Fig.6 we show the effect of changes in $`Q_o`$ and in the spectral index of the electrons in the input events, $`p^{}`$, on the overall electron population. The inputs for the dashed line are from Model 3. Line A shows the results of decreasing the spectral index using inputs from Model 6. Line B shows the results of decreasing $`Q_o`$ using inputs of Model 7. The results in this figure show that allowed variations in $`Q_o`$ and $`p^{}`$ are quite small since even modest changes in one of these parameters alone will result in an electron population that cannot produce the radio data.

We have investigated the extent to which we can mimic the best fit electron population shown in Model 3 by models which change $`p^{}`$ and $`Q_o`$ (the continuous injections). Inputs are shown in Model 8 and 9. The resultant electron distributions are shown in Fig.7. As can be seen there is only a limited amount of freedom in constructing this type of model since large changes in these offsetting parameters which can be used to fit the EUV data will quickly destroy the fit to the radio data. Good fits will be obtained if $`p^{}`$ is $``$ 2.32.

In Fig.8 we show the effect of different thermal electron densities on the cosmic ray electron population. The dashed line is the product of the inputs in Model 3; the solid line is the result using the inputs in Model 10. The electron densities in this model are an order of magnitude smaller than those in Model 3. As can be seen, the cosmic ray electron population is not greatly sensitive to changes in the thermal electron density in the higher energy range.

In Fig.9 we show the constraints imposed on the time of the initial event in our composite model. The input parameters are from Model 3. The dashed line corresponds to an evolutionary time of 1 Gyr, and the solid line is the electron population after 1.4 Gyr.

We also investigated the effect of a 5 $`\mu `$G field with regards to the hypothesis that both the EUV excess and the high energy X-ray flux are produced by the ICS mechanism with the same population of electrons. In this work we used Models 3, 8, 9, 10, and 11 which are capable of reproducing both the EUV and synchrotron radio flux. We investigated whether any of these models were also capable of producing the high energy X-ray flux by the ICS mechanism. The results are shown in Fig.10. In this figure, the top line is the observed excess hard energy flux (Rephaeli et al. 1999). The lower lines are the fluxes produced by the ICS mechanism in all five models. The flux produced by the ICS mechanism in all of these models is more than two orders of magnitude below the observed flux. Changes in the input characteristics of the injection events in these different models do not affect the predicted flux significantly.

## 4 DISCUSSION AND CONCLUSIONS

Atoyan & Völk (2000) are the only authors who have explored the production of an EUV excess in a high ($`1\mu `$G) intracluster field. They explored a model which is different in formal concept from those considered here, but is similar in form to our model with a large primary event followed by smaller injections. In specific, they considered two populations of cosmic rays. One is an old, highly evolved population which has undergone re-acceleration and compression; this population produces the EUV emission. A second higher energy population of cosmic rays which has been produced more recently by secondary origin is invoked to produce the radio emission. While this model may be appropriate, we note that models invoking secondary electrons for the radio emission have not been successful in producing the pronounced steeping of the spectral index of the radio emission with increasing distance from the cluster center that has been observed in well-studied radio halos such as Coma and Perseus. It has not been demonstrated that the Atoyan & Völk model can surmount this obstacle.

In this work we have studied general classes of models to establish which of these are capable of producing both the EUV and the synchrotron radio emission in a $`5\mu `$G cluster field. We studied models with a single initial event, models with continuous injection events, and combinations of these two. We find that in all models with just a single initial event, either the EUV emission will be underproduced or the synchrotron radio emission will be overproduced depending on the size of the event. With the passage of time, synchrotron losses (which are dominant in the high field case) will remove the higher energy components of the initial electron population, but the distribution required to produce the synchrotron radio emission is not produced. This problem is basic to all models with a single electron injection event. Continuous models encounter a different problem. An electron spectrum can be constructed which will produce the synchrotron radio emission. However, this population is not capable of producing the EUV flux.

The combination of a primary event followed by smaller continuous events can produce the required electron distribution. However, the initial event must have occurred at a time no earlier than 1.4 Gyr in the past or the higher energy electrons will have been lost and insufficient EUV flux will be produced. The electron energy loss rate is sufficiently steep that even substantially increasing the number of electrons in the initial event will not substantially change this time constraint. On the other hand, the time of the initial event must be no more recent than 1 Gyr ago or the synchrotron radio emission will be overproduced. This requires that the injection event must have occurred in the time span 1.4 Gyr $`<`$ injection $`<`$ 1 Gyr.

Finally, we investigated the effects of a 5$`\mu `$G field if both the EUV excess and the high energy X-ray flux are produced by ICS mechanism with the same population of electrons. The results shown in Figure 10 show that the high energy X-ray emission produced via the ICS process falls far short of the observed flux for all of these models. We conclude it is impossible to produce both the EUV flux and the hard X-ray emission by ICS from the same population of electrons with reasonable electron injection scenarios. We note that Atoyan & Völk (2000) have reached similar conclusions.

The project is supported by the National Science Council of Taiwan (grant number NSC 90-2112-M-008-040) and the Graduate Institute of Astronomy, NCU. S.B. thanks the faculty and staff of the Graduate Institute of Astronomy, NCU, and the National Center for Theoretical Science of Taiwan for their support and warm hospitality during an extended stay.


text:
# A new interpretation of the QCD phase transition and of strangeness as QGP signature

## 1 Introduction

It has been repeateadly demonstrated in the literature  that in many cases the final state of nuclear and particle collisions is compatible with the hypothesis of an equilibrated source. In the following the discussion will concern only colliding systems for which this finding holds (measured by the $`\chi ^2`$ of thermal model fits). Much work is presently concentrated on identifying the QCD phase transition at a certain $`\sqrt{s}`$, separating the colliding systems which go through the phase transition from those which dont. This defines the problem that we address. One way to approach this problem is to compare the estimated thermodynamic parameters temperature ($`T`$) and baryon chemical potential ($`\mu _B`$) at the chemical freeze-out extracted by thermal models from the data, to the expected critical ($`T_c`$, $`\mu _B(c)`$) values of QCD (e.g. ). However, the theoretical expectations for the critical temperature $`T_c`$ are uncertain. It is thus obviously interesting to try to extract informations on the critical parameters from the data, without any use of the QCD predictions. In this spirit, we approach this problem in the present study in a different way, namely by extrapolating the data to equivalent systems with zero $`\mu _B`$ and studying their $`T`$ . To see why this is interesting, we follow a line of arguments : If all measured colliding systems above a certain $`\sqrt{s}`$ are heated enough to go through the phase transition and back, they would appear to have the same $`T`$ (if they don’t freeze-out in a considerably different way) which we note as ’limiting’ hadronic temperature $`T_{lim}`$. The colliding systems which do not reach the $`T_c`$ at any time because of their small $`\sqrt{s}`$, will exhibit a final $`T`$ smaller than $`T_{lim}`$. Therefore, in order to separate the colliding systems which went through the transition from the ones which did not, one can investigate the $`T`$ of all colliding systems at conditions of $`\mu _B`$=0, as a function of $`\sqrt{s}`$ or $`dN/dy`$ at midrapidity, and search for an increase of $`T`$ followed by a saturation starting at the $`\sqrt{s}`$ where the $`T_c`$ is reached. This would work if we use always the same projectile and target combination. Otherwise, one should correct for the fact that different projectile/target combinations at the same $`\sqrt{s}`$, reach different initial energy densities. For this reason, we will investigate the $`T`$ as a function of the initial energy density ($`ϵ_i`$) reached at each collision after 1 fm/c based on the Bjorken estimate  and also using other Ansätze especially at low $`\sqrt{s}`$ . For a discussion on the uncertainty on $`ϵ_i`$ see . The same conclusions can however be reached by the present analysis, while looking only at one projectile/target combination A+A, namely with $`A200`$ (that is at Pb+Pb and Au+Au central collisions), as a function of e.g. $`\sqrt{s}`$ instead of $`ϵ_i`$.

## 2 Energy density dependence of temperature and $`\lambda _s`$

1. $`T`$: We compare the ratios of experimentally measured hadron yields in nuclear collisions with the prediction of a thermal model of non interacting free hadron resonances (for details see ). We extract the thermodynamic parameters describing best the particle source: temperature, baryochemical potential ($`\mu _B`$) and strangeness chemical potential ($`\mu _S`$), imposing exact strangeness conservation. We then extrapolate all thermodynamic states with parameters ($`T`$, $`\mu _B`$, $`\mu _s`$) to equivalent states at zero chemical potentials ($`T(\mu _B=0)`$, 0, 0) along isentropic paths. As demonstrated in figure 1 the resulting temperature at $`\mu _B=0`$ rises and saturates above $`ϵ_i`$ 1 $`GeV/fm^3`$. We interpret the onset of saturation of $`T`$ as due to the reach of an initial temperature greater or equal to $`T_c`$. We therefore extract from the onset of saturation in fig. 1 the ’critical’ $`\sqrt{s}`$ of $``$ 8.8 GeV/nucleon pair respectively the critical $`ϵ_c`$ of $``$ 1 +0.3 -0.5 GeV/fm3 of the QCD phase transition, independently of the QCD predictions. An important assumption needed to support this interpretation is that the cooling until freeze-out of the particle source, is not significantly (as compared to the errors) dependent on the $`\sqrt{s}`$, in the range between $`\sqrt{s}`$=2 and 9 GeV for A+A collisions. 111Note that a similar behaviour (rise and saturation) was found in the $`ϵ_i`$ dependence of the kaon number density $`\rho _k`$ , which may be related to the behaviour seen in fig. 1. A discussion of $`T`$ rise and saturation can be found in , but not for the $`\mu _B`$=0 case. However, it is only when investigating systems with the same $`\mu _B`$, where this behaviour is expected to be exact. The importance of using a common $`\mu _B`$ , becomes more apparent when investigating strangeness. This is the next topic.

2. $`\lambda _s`$: The so called ’strangeness suppression factor’ $`\lambda _s=\frac{(2\overline{s})}{(\overline{u}+\overline{d})}`$ at nonzero $`\mu _B`$, rises with $`ϵ_i`$ up to 1 GeV/fm3 (line (a) in fig. 2, left) and then decreases (line (b) in fig. 2, left). After extrapolating to $`\mu _B=0`$, $`\lambda _s`$ rises and saturates universally above $`ϵ_i`$ $``$ 1 GeV/fm3 (fig. 2, right). Therefore the peak of $`\lambda _s`$ at the 40 A GeV Pb+Pb (point at $`ϵ_i`$ = 1 GeV/fm3 in fig. 2, left) is due to the high $`\mu _B`$ reached there. Furthermore, the so called ’strangeness enhancement’ defined usually as double ratio of strange/pion ratio in A+A over p+p collisions at the same $`\sqrt{s}`$, can be traced back to the different $`\mu _B`$ and $`ϵ_i`$ of those reactions. In particular this explains why the so defined strangeness enhancement increases with decreasing $`\sqrt{s}`$ (fig. 7 in ), since this corresponds to increasing $`\mu _B`$. If one eliminates the bias due to different $`\mu _B`$ by considering only $`\lambda _s`$ at $`\mu _B=0`$, the non-trivial -with respect to the phase transition- ’strangeness enhancement’ is revealed (fig. 2, right) as a consequence of the behaviour of $`T`$ seen in fig. 1. We therefore redefine the ’strangeness enhancement’, as an enhancement of $`\lambda _s`$ in all thermalised systems above $`ϵ_c`$, as compared to all thermalised systems with $`ϵ_i<ϵ_c`$, at the same $`\mu _B`$. 222Discussions of the QCD phase transition appearing possibly between AGS and SPS energy can be found e.g. in , however not at a common $`\mu _B`$. We also analysed for this conference ratios from $`e^+e^{}`$ collisions at $`\sqrt{s}`$=3.6 GeV  and found a T of 124 +15 -20 MeV, and $`\lambda _s`$=0.25 +0.06 -0.09 with $`\chi ^2/DOF`$=4 10-3/1. This T is below $`T_{lim}`$, however the $`ϵ_i`$ is difficult to extract due to lack of data. Extrapolating the A+A data in fig. 2, left, to higher $`ϵ_i`$, we find that the approximately netbaryonfree limit of $`\lambda _s`$ in Au+Au collisions, is expected to be reached at LHC energies .

## Conclusions

The problem we address is how to identify the onset of the QCD phase transition using measured light (u,d,s)-flavoured hadrons, and separate the colliding systems which go through the QCD phase transition from the ones which dont, without any use of the QCD predictions, as well as the role of strangeness. We achieve this by estimating for the first time  the equivalent temperature at zero baryochemical potential of the thermodynamic states describing several measured particle and nuclear collisions. We find at zero chemical potential ($`\mu _B`$=0) a universal rise and saturation of both the $`T`$ and of the strangeness suppression factor $`\lambda _s`$ (=$`\frac{2\overline{s}}{\overline{u}+\overline{d}}`$) with increasing initial energy density ($`ϵ_i`$) (fig. 1 and fig. 2, right). The onset of saturation of both $`T`$ and $`\lambda _s`$ at $`\mu _B=0`$, allows to discriminate systems which go through the QCD phase transition from those which dont. The critical energy density is thus estimated at the onset of saturation to be $`ϵ_c`$ $``$ 1 +0.3 -0.5 GeV/fm3, corresponding approximately to a $`\sqrt{s}`$ of $``$ 8.8 GeV for central Pb+Pb collisions. Further, we identify trivial and non-trivial (with respect to the phase transition) sources of ’strangeness enhancement’: e.g. the maximum of $`\lambda _s`$ at the 40 A GeV Pb+Pb collisions , and other phenomena of ’strangeness enhancement’, are ’trivially’ traced back to the different baryochemical potentials and $`ϵ_i`$ of the compared systems (compare fig. 2 right with $`\mu _B`$=0 and left with nonzero $`\mu _B`$). This explaines why e.g. the strange/pion double ratios in A+A over p+p collisions both at the same $`\sqrt{s}`$, increase with decreasing $`\sqrt{s}`$, because $`\mu _B`$ increases and $`ϵ_i`$ differs for A+A and p+p collisions at the same $`\sqrt{s}`$. It could explain also an enhancement of $`\lambda _s`$ in p+A over p+p collisions. We redefine the non-trivial ’strangeness enhancement’, as the increase of $`\lambda _s`$ in all (thermalised) systems which reached $`ϵ_c`$ as compared to all (thermalised) systems which did not, at the same $`\mu _B`$. It is a consequence of the $`ϵ_i`$ dependence of $`T`$. The netbaryonfree $`\lambda _s`$ limit is estimated to be approximately reached in Au+Au collisions at the LHC . In conclusion, in contrast to external non-equilibrium signatures like charmonia suppression or jet quenching which are important but may be overcritical, hadrons with $`u,d,s,\overline{u},\overline{d},\overline{s}`$ quarks play a special role as equilibrium signatures and part of the plasma itself, in allowing to extract the critical parameters of the QCD phase transition.

## Acknowledgments

We thank the Schweizerischer Nationalfonds for their support, as well as K. Pretzl and P. Minkowski for fruitfull discussions. We also thank the organisers of ISMD2001 for creating an open and scientifically fruitfull atmosphere.


text:
# Spin polarized neutron matter and magnetic susceptibility within the Brueckner–Hartree–Fock approximation

## I Introduction

The study of the magnetic properties of dense matter is of considerable interest in connection with the physics of pulsars. These objects, since the suggestion of Gold , are generally believed to be rapidly rotating neutron stars with strong surface magnetic fields of the order of $`10^{12}`$ Gauss. Several authors have studied the possible existence of a phase transition to a ferromagnetic state on pure neutron systems at densities corresponding to the theoretically stable neutron stars. Brownell and Callaway , and Rice  considered a hard sphere gas model and showed that the ground state of the neutron gas becomes ferromagnetic at $`k_F2.3`$ fm-1. Silverstein  and Østgaard  found that the inclusion of long range attraction significantly increased the ferromagnetic transition density (e.g. Østgaard predicted the transition to occur at $`k_F4.1`$ fm-1 using a simple central potential with hard core only for singlet spin states). Clark  and Pearson and Saunier  calculated the magnetic susceptibility for low densities ($`k_F2`$ fm-1) using more realistic interactions. Pandharipande et al. , using the Reid soft-core potential, performed a variational calculation arriving to the conclusion that such a transition was not to be expected for $`k_F5`$ fm-1. Early calculations of the magnetic susceptibility within the Brueckner theory were performed by Bäckmann and Källman  employing the Reid soft–core potential, and results from a Correlated Basis Function calculation were obtained by Jackson et al.  with the Reid $`v_6`$ interaction. A different point of view was followed by Vidaurre et al. , who employed neutron-neutron effective interactions of Skyrme type, finding the ferromagnetic transition at $`k_F1.73`$$`1.97`$ fm-1.

In connection with the problem of the neutrino diffusion in dense matter, Fantoni et al.  have recently employed a new quantum simulation technique (the so-called Auxiliary Field Diffusion Monte Carlo method (AFDMC)) using realistic interactions (based upon the Argonne $`v_{18}`$ two–body potential  plus Urbana IX three–body potential ), and have found that the magnetic susceptibility of neutron matter shows a strong reduction of about a factor $`3`$ with respect its Fermi gas value. They pointed out that such a reduction may have strong effects on the mean free path of a neutrino in dense matter and, therefore, it should be taken into account in the studies of supernovae and proto-neutron stars.

In this work we employ the Brueckner–Hartree–Fock (BHF) approximation, using the realistic Nijmegen II and Reid93  nucleon-nucleon interacions, to study spin polarized neutron matter properties such as the total energy per particle and the magnetic susceptibility. We employ the so-called continuous prescription when solving the Bethe–Goldstone equation. As shown by Song et al. , the effects from three-body clusters are dimished in this prescription. We also explore in this work the dependence of the total energy per particle on the spin polarization, finding that up to a very good approximation this dependence is parabolic.

The paper is organized in the following way. In Section II the theoretical background of our calculation is briefly reviewed. The construction of the neutron-neutron $`G`$-matrices and the calculation of the total energy per particle are shown in Section II A, whereas the magnetic susceptibility is determined in Section II B. Our results are presented in Section III. Finally, a short summary and the main conclusions are given in Section IV.

## II Theoretical background

In this section we briefly show how to evaluate, in the BHF approximation, the total energy per particle and the magnetic susceptibility of a system of neutrons in which we assume that the density of particles with spin up, $`\rho ^{}`$, is different from that with spin down, $`\rho ^{}`$.

### A Energy per particle

Our calculation of the total energy per particle starts with the construction of $`G`$-matrices, which describe in an effective way the interaction between two neutrons, for each one of the spin combinations ($``$, $``$, $``$ or $``$), in the presence of a surrounding medium. They can be obtained by solving the following integral Bethe–Goldstone equations

$`G_,=V_,+V_,{\displaystyle \frac{Q_,}{\omega ϵ_{}ϵ_{}+i\eta }}G_,,`$ (1)
$`G_,=V_,+V_,{\displaystyle \frac{Q_,}{\omega ϵ_{}ϵ_{}+i\eta }}G_,,`$ (2)
$`\left(\begin{array}{cc}G_,& G_,\\ G_,& G_,\end{array}\right)=\left(\begin{array}{cc}V_,& V_,\\ V_,& V_,\end{array}\right)+(\begin{array}{cc}V_,& V_,\\ V_,& V_,\end{array})(\begin{array}{cc}\frac{Q_,}{\omega ϵ_{}ϵ_{}+i\eta }& 0\\ 0& \frac{Q_,}{\omega ϵ_{}ϵ_{}+i\eta }\end{array})\left(\begin{array}{cc}G_,& G_,\\ G_,& G_,\end{array}\right)`$

In the above expressions the first (last) two subindices indicate the spin projections of the two neutrons in the initial (final) state, $`V`$ is the bare nucleon-nucleon interaction, $`Q`$ is the Pauli operator which allows only intermediate states compatible with the Pauli principle, and $`\omega `$ is the starting energy defined as the sum of single–particle energies, $`ϵ_{()}`$, of the interacting neutrons. Note that $`G_,`$ and $`G_,`$ are obtained from a coupled channel equation due to the mixing induced by the interaction. One can equivalently solve the Bethe–Goldstone equation in the spin–coupled basis, where the interaction is diagonal, although in that case the Pauli operator is non-diagonal. However, the conventional angle-average of the Pauli operator makes it diagonal, thus reducing the problem to an uncoupled one in each total spin channel.

The single–particle energy of a neutron with momentum $`k`$ and spin projection $`\sigma =`$($``$) is given by

$$ϵ_\sigma =\frac{\mathrm{}^2k^2}{2m}+U_\sigma (k),$$
(3)
where the single–particle potential $`U_\sigma (k)`$ represents the average field felt by the neutron due to its interaction with the other neutrons of the system. In the BHF approximation it is given by

$$U_\sigma (k)=\text{Re}\underset{\sigma ^{}=,}{}\underset{k^{}k_F^\sigma ^{}}{}\stackrel{}{k}\stackrel{}{k^{}}|G_{\sigma \sigma ^{},\sigma \sigma ^{}}(\omega =ϵ_\sigma +ϵ_\sigma ^{})|\stackrel{}{k}\stackrel{}{k}^{}_𝒜,$$
(4)
where a sum over the two Fermi seas of spin up and down, characterized by $`k_F^{}=(6\pi ^2\rho ^{})^{1/3}`$ and $`k_F^{}=(6\pi ^2\rho ^{})^{1/3}`$ respectively, is performed and the matrix elements are properly antisymmetrized.

Once a self-consistent solution of Eqs. (2) and (4) is obtained, the total energy per particle is easily calculated

$$\frac{E}{N}=\underset{\sigma =,}{}\underset{kk_F^\sigma }{}\left(\frac{\mathrm{}^2k^2}{2m}+\frac{1}{2}U_\sigma (k)\right).$$
(5)
This quantity is a function of $`\rho ^{}`$ and $`\rho ^{}`$ or, equivalently, of the total density $`\rho =\rho ^{}+\rho ^{}`$ and the spin polarization $`\mathrm{\Delta }`$, defined as

$$\mathrm{\Delta }=\frac{\rho ^{}\rho ^{}}{\rho }.$$
(6)
Note that the value $`\mathrm{\Delta }=0`$ corresponds to non-polarized or paramagnetic ($`\rho ^{}=\rho ^{}`$) neutron matter, whereas $`\mathrm{\Delta }=\pm 1`$ means that the system is totally polarized, i.e., all the spins are aligned in the same direction.

### B Magnetic susceptibility

The magnetic susceptibility of a system characterizes the response of this system to a magnetic field and gives a measure of the energy required to produce a net spin alignment in the direction of the field. It is defined as

$$\chi =\left(\frac{}{}\right)_{=0},$$
(7)
where $``$ is the magnetization of the system per unit volume given by

$$=\mu (\rho ^{}\rho ^{})=\mu \rho \mathrm{\Delta },$$
(8)
with $`\mu `$ the magnetic moment of a neutron, and $``$ is the magnetic field which can be obtained from

$$=\rho \left(\frac{(E/N)}{}\right)_{=0}=\frac{1}{\mu }\left(\frac{(E/N)}{\mathrm{\Delta }}\right)_{\mathrm{\Delta }=0}.$$
(9)
Using Eqs. (8) and (9), the magnetic susceptibility can be written as

$$\chi =\frac{\mu ^2\rho }{\left(\frac{^2(E/N)}{\mathrm{\Delta }^2}\right)_{\mathrm{\Delta }=0}},$$
(10)
where the second derivative can be taken at $`\mathrm{\Delta }=0`$ if the field $``$ is assumed to be small.

It is customary to study the magnetic susceptibility in terms of the ratio $`\chi /\chi _F`$, where $`\chi _F`$ is the magnetic susceptibility of a free Fermi gas, usually known as Pauli susceptibility. It can be straightforwardly obtained from Eq. (10) and the total energy per particle of the free Fermi gas

$$\chi _F=\frac{\mu ^2m}{\mathrm{}^2\pi ^2}k_F.$$
(11)
where the Fermi momentum $`k_F=(3\pi ^2\rho )^{1/3}`$ is related to $`k_F^{}`$ and $`k_F^{}`$ through the relations

$$\begin{array}{c}k_F^{}=k_F(1+\mathrm{\Delta })^{1/3}\\ k_F^{}=k_F(1\mathrm{\Delta })^{1/3}.\end{array}$$
(12)
## III Results

The total energy per particle for totally polarized (solid lines) and non-polarized (dashed lines) neutron matter is shown in Fig. 1 as a function of the density. Results for the Nijmegen II interaction are plotted on the left panel, whereas those corresponding to the Reid93 interaction are shown on the right panel. As can be seen from the figure, for both interaction models, the energy of totally polarized neutron matter is always more repulsive than non-polarized neutron matter in all the density range explored. This additional repulsion can be understood, firstly, in terms of the kinetic energy contribution, which is larger in the totally polarized case than in the non-polarized one; and secondly, in terms of the potential energy contribution because, due to symmetry arguments, all partial waves with even orbital angular momentum $`L`$ (some of them attractive, as the important $`{}_{}{}^{1}S_{0}^{}`$) are excluded in totally polarized neutron matter. In order to illustrate this, we have plotted in Fig. 2 the separate kinetic (left panel) and potential (right panel) energy contributions for the Nijmegen II interaction model (similar results are obtained for the Reid93 one, but they are not included in order to make the discussion more clear). An interesting conclusion which can be inferred from these results is that a phase transition to a ferromagnetic state is not to be expected from our calculation. If such a transition would exist a crossing of the energies of the totally polarized and the non-polarized cases would be observed at some density, indicating that the ground state of the system would be ferromagnetic from that density on. As can be seen in Fig. 1, there is no sign of such a crossing and, on the contrary, it becomes less favourably as the density increases.

We have shown results for totally polarized and non-polarized neutron matter. Let us consider now an intermediate situation in which not all the spins, but a part of them, are aligned in a given direction, and let us examine the dependence of the total energy per particle in the spin polarization $`\mathrm{\Delta }`$. This dependence is shown in Fig. 3 for five different densities ($`\rho _0/2,\rho _0,2\rho _0,5\rho _0`$ and $`7\rho _0`$, being $`\rho _0=0.17`$ fm-3 the saturation density of nuclear matter). As in Fig. 1, results for the Nijmegen (Reid93) are shown on the left (right) panel. Circles, squares, diamonds and triangles correspond to our BHF results, whereas solid lines correspond to the parabolic approximation discussed below. As can be seen from this figure, and as it was expected, $`E/N`$ is symmetric in $`\mathrm{\Delta }`$. It can also be seen in this figure that $`E/N`$ shows a minimum at $`\mathrm{\Delta }=0`$ for all the densities considered, being this again an indication that the ground state of neutron matter is paramagnetic. Another interesting thing is to note that this dependence is up to a very good approximation parabolic, being this parabolic character only slightly lost at very large densities. Therefore, in the same spirit as it is done in nuclear matter to determine the symmetry energy, one can try to characterize this dependence in the following simple analytic form

$$\frac{E}{N}(\rho ,\mathrm{\Delta })=\frac{E}{N}(\rho ,0)+a(\rho )\mathrm{\Delta }^2,$$
(13)
where, assuming the quadratic dependence to be valid up to $`\mathrm{\Delta }=1`$ as our results indicate, the value of $`a(\rho )`$ can be easily obtained for each density as the difference between the total energy per particle of totally polarized and non-polarized neutron matter

$$a(\rho )=\frac{E}{N}(\rho ,\pm 1)\frac{E}{N}(\rho ,0).$$
(14)
The magnetic susceptibility can be evaluated in a very simple way if the parabolic dependence of Eq. (13) is assumed, giving

$$\chi (\rho )=\frac{\mu ^2\rho }{2a(\rho )}.$$
(15)
In Fig. 4 the ratio $`\chi /\chi _F`$ is shown as a function of the density. The solid line shows the result for the Nijmegen II interaction, whereas the dashed line corresponds to the one obtained with Reid93. Starting from 1, the ratio decreases rapidly for small densities and more slowly as density increases. It can be inferred again from this figure that a ferromagnetic phase transition, which would be signaled by an infinite discontinuity giving rise to a change of sign in $`\chi /\chi _F`$, is not seen and not expected at larger densities either.

Finally, our results for $`\chi /\chi _F`$ are compared in Table I with those of the recent calculation performed by Fantoni et al. , shown in columns labelled AU6’ and AU8’. As can be seen from the table, there is a very good agreement between these results and ours. For completeness, we show in parentheses the results obtained when the standard discontinuous prescription is used in solving the Bethe–Goldstone equation. In both prescriptions, the results are very similar, which is not surprising due to the fact that $`\chi `$ is obtained from an energy difference (see Eqs. (14) and(15)) which partly cancels the possible discrepancies. Only for densities larger than $`2\rho _0`$, the discontinuous prescription results differ more than $`10\%`$ from the continuous ones.

## IV Summary and Conclusions

Employing realistic modern nucleon-nucleon interactions (Nijmegen II and Reid93) we have performed a Brueckner–Hartree–Fock calculation of spin polarized neutron matter properties. We have studied the total energy per particle of neutron matter as a function of the density and the spin polarization $`\mathrm{\Delta }`$. We have found that in the range of densities explored (up to $`7\rho _0`$) totally polarized matter is always more repulsive than non-polarized matter, being this an indication that a phase transition of the system to a ferromagnetic state is not expected.

We have seen that the total energy per particle is not only symmetric on the spin polarization $`\mathrm{\Delta }`$, as it was expected, but also parabolic in a very good approximation up to $`\mathrm{\Delta }=1`$ even at high densities. This finding supports the calculation of the magnetic susceptibility of neutron matter by using only the energies of the spin symmetric and fully polarized systems.

Finally, we have calculated the magnetic susceptibility of the system as a function of the density, finding a very good agreement with a recent Monte Carlo calculation .

## Acknowledgements

This work is partially supported by the DGICYT contract No. PB98-1247 and by the Generalitat de Catalunya grant No. 2000SGR00024 (Spain). The authors acknowledge useful discussions with Prof. J. Navarro. One of the authors (I.V.) wishes to acknowledge support from a post-doctoral fellowship of the Istituto Nazionale di Fisica Nucleare (Italy).