text:
# 1 Introduction

## 1 Introduction

The instability of planar films against coherently strained three-dimensional (3D) islands in highly mismatched epitaxy is a subject of intense research in recent time owing to their possible optoelectronic applications as quantum dots. The term “coherent Stranski-Krastanov (SK) growth” has been coined for this case of formation of 3D islands that are strained to fit the underlying wetting layer at the interface but are largely strain-free near their top and side walls. This term was introduced in order to distinguish this case from the “classical” SK growth in which the lattice misfit is accommodated by misfit dislocations at the interface.

Experimental studies of arrays of coherent 3D islands in SK growth of highly mismatched semiconductor materials one on top of the other have shown surprisingly narrow size distributions of the islands, (see also Ref.  and the references therein). (It is worth noting that a narrow size distribution has been established also in the Volmer-Weber growth of metals on insulators in the absence of a wetting layer.) This phenomenon, known in the literature as self-assembly (for a review see Ref. ), is highly desirable as it guarantees a specific optical wavelength of the array of quantum dots. The physics of this self-assembly is still not understood in spite of the numerous thermodynamic and kinetic studies. For example, Priester and Lannoo found that two-dimensional (2D) islands with a monolayer height act as precursors of the 3D pyramidal islands, (see also Ref. ). The energy per atom of the 2D islands possesses a minimum for a certain volume, but the 3D islands become energetically favored at a smaller size. Thus, at some critical surface coverage, the 2D islands spontaneously transform into 3D islands preserving a nearly constant volume during the 2D-3D transformation. The resulting size distribution reflects that of the 2D islands which is very narrow. This picture has been recently corroborated by Ebiko et al. who found that the volume distribution of InAs/GaAs self-assembled quantum dots agrees well with the scaling function that is characteristic for the two-dimensional submonolayer homoepitaxy.  Korutcheva et al. and Markov and Prieto  reached the same conclusion with the exception that the 2D-3D transformation was found to take place through a series of intermediate states with discretely increasing thickness (one, two, three, etc. monolayers-thick islands) that are stable in separate consecutive intervals of volume. Khor and Das Sarma arrived to the same conclusion by using Monte Carlo simulations.

In a recent paper, Prieto and Markov discussed the formation of coherent 3D islands within the framework of the traditional concept of wetting. As is well known the wetting parameter which accounts for the energetic influence of a crystal B in the heteroepitaxial growth of a crystal A on top of it is defined as (for a review see Ref. )

$$\mathrm{\Phi }=1\frac{E_{\mathrm{AB}}}{E_{\mathrm{AA}}}$$
(1)
where $`E_{\mathrm{AA}}`$ and $`E_{\mathrm{AB}}`$ are the energies per atom required to disjoin a half-crystal A from a like half-crystal A and from an unlike half-crystal B, respectively.

The mode of growth of a thin film is determined by the difference $`\mathrm{\Delta }\mu =\mu (n)\mu _{3\mathrm{D}}^0`$, where $`\mu (n)`$ and $`\mu _{3\mathrm{D}}^0`$ are the chemical potentials of the film (as a function of its thickness $`n`$) and of the bulk material A, respectively. The chemical potential of the bulk crystal A is given at zero temperature by the work $`\varphi _{\mathrm{AA}}`$ (taken with a negative sign) to detach an atom from the well known kink or half-crystal position. The latter name is due to the fact that an atom at this position is bound to a half-atomic row, a half-crystal plane and a half-crystal block. In the case of a monolayer-thick film of A on the surface of B the chemical potential of A is given by the analogous work $`\varphi _{\mathrm{AB}}`$ with the exception that the underlying half-crystal block of A is replaced by a half-crystal block of B. Thus $`\mathrm{\Delta }\mu =\varphi _{\mathrm{AA}}\varphi _{\mathrm{AB}}`$. In the simplest case of additivity of bond energies the difference $`\varphi _{\mathrm{AA}}\varphi _{\mathrm{AB}}`$ reduces to $`E_{\mathrm{AA}}E_{\mathrm{AB}}`$ as the lateral bondings cancel each other. Then $`\mathrm{\Delta }\mu `$ is proportional to $`\mathrm{\Phi }`$, i.e. $`\mathrm{\Delta }\mu =E_{\mathrm{AA}}\mathrm{\Phi }`$. It follows that it is the wetting parameter $`\mathrm{\Phi }`$ which determines the mechanism of growth of A on B. In the two limiting cases of growth of isolated 3D islands of A directly on the surface of B \[Volmer-Weber growth, characterized by incomplete wetting ($`0<\mathrm{\Phi }<1`$) and any misfit $`\epsilon _0=\mathrm{\Delta }a/a`$\], or by consecutive formation of monolayers of A on B \[Frank - van der Merwe growth, with complete wetting ($`\mathrm{\Phi }0`$) and $`\epsilon _00`$\], $`\mathrm{\Delta }\mu `$ goes asymptotically to zero from above and from below, respectively, but it changes sign in the case of growth of 3D islands of A on a thin wetting layer of A on the substrate B \[Stranski-Krastanov growth, complete wetting at the beginning ($`\mathrm{\Phi }0`$) and $`\epsilon _00`$\]. Equation $`\mathrm{\Delta }\mu =E_{\mathrm{AA}}\mathrm{\Phi }`$ is thus equivalent to the familiar 3-$`\sigma `$ criterion of Bauer.

The Stranski-Krastanov morphology appears as a results of the interplay of the film-substrate bonding, misfit strain and the surface energies. A wetting layer with a thickness of the order of the range of the interatomic forces is first formed (owing to the interplay of the A-B interaction and the strain energy accumulation) on top of which partially or completely relaxed 3D islands nucleate and grow. The 3D islands and the thermodynamically stable wetting layer represent necessarily different phases. If this were not the case, the growth would continue by 2D layers. Thus we can consider as a useful approximation to regard the 3D islanding on top of the uniformly strained wetting layer as a Volmer-Weber growth. That requires the mean adhesion of the atoms that belong to the base plane of the 3D islands to the stable wetting layer to be smaller than the cohesion between them. In other words, the wetting of the underlying wetting layer by the 3D islands must be incomplete. Otherwise, 3D islanding will not occur. In the Volmer-Weber growth the incomplete wetting is due mainly to the difference in bonding ($`E_{\mathrm{AB}}<E_{\mathrm{AA}}`$), the supplementary effect of the lattice misfit being usually smaller. In the coherent SK growth ($`E_{\mathrm{AB}}E_{\mathrm{AA}}`$), the incomplete wetting is due to the lattice misfit, which leads to the atoms at the edges of the islands to displace from the bottoms of the potential troughs provided by the atoms in the layer underneath.

As has been shown elsewhere, it is the incomplete wetting which determines the formation of dislocation-free 3D islands on top of the wetting layer in the case of coherent SK growth. In this case, however, the relation between $`\mathrm{\Delta }\mu `$ and $`\mathrm{\Phi }`$ is not as simple as given above. The bond energies are generally not additive, the misfit strain is relaxed mostly near the side and top walls and increasing the island’s thickness leads to larger displacements of the edge atoms from the bottoms of the potential troughs provided by the wetting layer and in turn to a decrease of the wetting. For this reason, in this work we define the wetting parameter $`\mathrm{\Phi }`$ as the difference of the interaction energies with the wetting layer of misfitting and non-misfitting 3D islands.

In the present paper we study the behavior of $`\mathrm{\Phi }`$ for islands which belong to an array of islands. We study the effect of the density of the array (the distance to the nearest neighbor islands), the size distribution (the difference in size of the neighboring islands), and the shape distribution (the slope of the side walls of the neighboring islands) on the wetting parameter $`\mathrm{\Phi }`$ of the considered island.

## 2 Model

We consider an atomistic model in $`1+1`$ dimensions (lateral size + height) which we treat as a cross section of the real $`2+1`$ dimensional case. An implicit assumption is that in the real $`2+1`$ dimensional model the monolayer islands have a compact rather than a fractal shape and that the lattice misfit is the same in both orthogonal directions. The 3D islands are represented by linear chains of atoms stacked one upon the other. Each upper chain is shorter than the lower one. The shape of the islands in our model is given by the slope of the side walls. For, example, an island which consists of consecutive chains with N, N-1, N-2… atoms has a 60 slope of the side walls, whereas an island with chains consisting of, say, N, N-5, N-10… atoms has a slope of 19.1, etc., where N is the number of atoms in the base chain. The array in the 1+1 dimensional space is represented by a row of 3 or 5 islands on a wetting layer consisting of several monolayers (Fig. 1). The distance between two neighboring islands is given by the number $`n`$ of vacant atomic positions betweeen the ends of their base chains and can be varied from one to infinity.

In order to simplify the computational procedure the “wetting layer” in our model is in fact composed of several monolayers of the true wetting layer consisting of atoms of the overlayer material A plus several monolayers of the unlike substrate material B. This composite wetting layer has the atom spacing of the substrate material B as in the real case, but for the sake of simplicity the atom bonding is that of the overlayer material A. We believe that the latter does not introduce a perceptible error as the energetic influence of the substrate B is screened by the true wetting layer A. We expect that this approximation underestimates to some extent the wetting parameter $`\mathrm{\Phi }`$ by making the composite substrate a bit softer than the real one (the A-A bonding is weaker than the B-B bonding, for a later review see Ref. ). We found that beyond 10 monolayers the studied parameter $`\mathrm{\Phi }`$ saturates its value. This is why in all cases given below, unless otherwise stated, we allowed 10 monolayers to relax.

As in our previous work (Ref. ), we make use of a simple minimization procedure. The atoms interact through a potential that can be easily generalized to vary its anharmonicity by adjusting two constants $`\mu `$ and $`\nu `$ ($`\mu >\nu `$) that govern separately the repulsive and attractive branches, respectively,

$`V(x)=V_o\left[{\displaystyle \frac{\nu }{\mu \nu }}e^{\mu (xb)}{\displaystyle \frac{\mu }{\mu \nu }}e^{\nu (xb)}\right],`$ (2)
where $`b`$ is the equilibrium atom separation. In this work, we have used $`\mu =2\nu `$ with $`\nu =6`$, which turns the potential (2) into the familiar Morse form.

Our program calculates the interaction energy of all the atoms as well as its gradient with respect to the atomic coordinates, i.e. the forces. Atoms in the islands and in the wetting layer are then allowed to relax. Relaxation of the system is performed by allowing the atoms to displace in the direction of the gradient in an iterative procedure until the forces fall below some negligible cutoff value. Periodic boundary conditions are applied in the lateral direction. We consider only interactions in the first coordination sphere in order to mimic the directional bonds that are characteristic for most semiconductor materials.

## 3 Results

Figure 2(a) shows the horizontal displacements of the atoms of the base chain from the bottoms of the potential troughs provided by the homogeneously strained wetting layer for a misfit of 7%. The considered island has two identical ones at a distance of $`n=5`$. This is the same behavior as predicted by the one-dimensional model of Frank and van der Merwe. The horizontal displacements increase with increasing island thickness (measured in number of monolayers) precisely as in the case of a rigid substrate and non-interacting islands. In contrast to the rigid substrate case, the vertical displacements of the edge atoms of the base chain of the islands and the underlying atoms of the uppermost monolayer of the wetting layer are directed downwards (Fig. 2(b)). It is worth noting that the same result has been found by Lysenko et al. in the case of homoepitaxial metal growth by using a computational method within the framework of the tight-binding model.

In spite of their downwards vertical displacements, the edge atoms are again more weakly bound to the underlying wetting layer, as in the rigid substrate models of Refs.  and  (Fig. 3). Shown in the same figure for comparison is also an island without neighbors. As seen, the edge atoms of the single island adhere more strongly to the substrate. This is in fact the essential physics behind the effect that the neighboring islands exert on the middle island. The central island looses to some degree contact with the substrate (in this case the wetting layer) and the wetting parameter is increased. We can interpret this as the wetting layer becoming stiffer under the influence of the neighboring islands.

The influence of the density of the array is demonstrated in Fig. 4. The values for 3 and 5 islands were calculated assuming equally spaced islands. These can be thus treated as a self-organized array. As expected the wetting parameter increases with decreasing distance between the islands or, in other words, with increasing array density.

Figure 5 shows the wetting parameter of the central island as a function of the size of the side islands. For this calculation, we considered three islands with the same thickness of 3 monolayers. Furthermore, both side islands have one and the same volume. Increasing the volume of the side islands leads to an increase of the elastic fields around them and to a further reduction of the bonding between the edge atoms of the central island and the wetting layer.

Figure 6 demonstrates one of the most important results, the effect of the size distribution on the wetting of the islands. It shows the behavior of the wetting parameter $`\mathrm{\Phi }`$ of the central island as a function of the number of atoms in the base chain (which is a measure of the volume) of the left island. In this case, the sum of the volumes (the total number of atoms) of the left and right islands is kept constant and precisely equal to the doubled volume of the central island. All three islands have the same thickness. The side walls of all the three islands make an angle of 60 with their bases. Thus the first point (and, by symmetry, also the last one) gives the maximum asymmetry in the size distribution of the array, the left island consists of 9 atoms whereas the right island is built of 105 atoms. The point at the maximum of wetting describes the monodisperse distribution – the three islands have one and the same volume of 57 atoms. As seen in the case of perfect self-assembly of the array the wetting parameter, or in other words, the tendency to clustering displays a maximum value.

The effect of the shape of the side islands, i.e. their facet angles, on the wetting parameter of the central island is demonstrated in Fig. 7. The slope of the facets of the central island is 60. The effect is greatest when the side islands have the steepest walls. The same result (not shown) is obtained when the central island has a different facet angle, e.g. 11. The explanation follows the same line as the one given above. The side islands with larger-angle side walls exert a greater elastic effect on the wetting layer and in turn on the displacements and the bonding of the edge atoms of the central island.

We studied also the stability of islands with a thickness increasing by one monolayer in the presence of two islands on both sides on a deformable substrate. The result in Fig. 8 shows the same behavior observed in Refs.  and  where rigid substrates were assumed. This means that the overall transformation from the precursor 2D islands to a macroscopic 3D islands takes place in consecutive stages in each of which the islands thicken by one monolayer. As shown in Refs.  and  the latter leads in turn to a critical misfit beyond which coherent 3D islanding takes place, and below which the lattice mismatch is accommodated by misfit dislocations. The existence of a critical misfit has been experimentally observed in a series of different systems. The energies computed in the case of the reference single islands always lie below the curves of the islands in an array. The difference obviously gives the energy of repulsion between the neighboring islands. It follows from the above that the presence of neighboring islands leads to a slight decrease of the critical misfit.

## 4 Discussion

For the discussion of the above results we have to bear in mind that a positive wetting parameter shows in fact a tendency of the deposit to form 3D clusters rather than a planar film. In the case of coherent SK growth, the non-zero wetting parameter is due to the weaker adhesion of the atoms that are closer to the islands edges. The presence of other islands, particularly with large angle facets, in the near vicinity of the considered island makes this effect stronger as seen in Fig. 3. The transformation of two-dimensional islands with a monolayer height into bilayer three-dimensional islands takes place by detachment of atoms from the edges and their subsequent jumping and collision on the top island’s surface. Thus this edge effect clearly demonstrates the influence of the lattice misfit on the rate of second layer nucleation and in turn on the kinetics of the 2D-3D transformation. The presence of neighboring islands facilitates and thus accelerates the formation of 3D clusters and their further growth. In a self-assembled population of islands the tendency to clustering is thus enhanced.

We can think of the flatter islands in our model (11 facet angle) as the famous “hut” clusters discovered by Mo et al., and of the clusters with 60 facet angles as the “dome” clusters. It is well known that clusters with steeper side walls relieve the strain much more efficiently than the flatter clusters (see the discussion in Ref. ; the planar film, which is the limiting case of the flatter islands with a facet angle equal to zero, does not relieve the strain at all.) We see that large-angle facet islands affect more strongly the growth of the neighboring islands, leading thus to a more narrow size distribution.

We further conclude that a self-assembled population of quantum dots is expected at comparatively low temperatures such that the critical wetting-layer thickness for 3D islanding to take place approaches an integer number of monolayers. In InAs/GaAs quantum dots, the reported values of the critical thickness were found to vary from 1.2 to 2 monolayers, (see also the discussion in Ref.  and the references therein). The critical wetting-layer thickness should be given by an integer number of monolayers plus the product of the 2D island density and the critical volume (or area) $`N_{12}`$. The 2D island density increases steeply with decreasing temperature. In such a case, a dense population of 2D islands will overcome simultaneously the critical size $`N_{12}`$ to produce 3D bilayer islands. The latter will interact maximally with each other from the very beginning of the 2D-3D transformation giving rise to a maximum wetting parameter and, in turn, to large-angle facets and a narrow size distribution. This is in agreement with the observations of Le Tanh et al. in the case of nucleation and growth of self-assembled Ge quantum dots on Si(001). At 700C, a population of islands with a concentration of the order of 1$`\times 10^7`$ \- 1$`\times 10^8`$ cm<sup>-2</sup> is obtained; the islands have the shape of a truncated square pyramid with four side wall facets formed by (105) planes with an inclination angle of about 11 and the size distribution of the islands is quite broad. On the other hand, at 550C, a population of islands with an areal density of the order of 1$`\times 10^9`$ \- 1$`\times 10^{10}`$ cm<sup>-2</sup> is observed, the islands have larger angle (113) facets and their size distribution is much more narrow.

In summary, we have shown that the presence of neighboring islands decreases the wetting of the substrate (in this case the wetting layer) by the 3D islands. The larger the density of the array, the weaker the wetting. Neighboring islands with steeper side walls reduce more strongly the wetting of the considered island. The wetting parameter displays a maximum (implying a minimal wetting) when the array shows a monodisperse size distribution. We should expect optimum self-assembled islanding at lower temperatures such that the 2D-3D transformation takes place at the maximum possible island density.

###### Acknowledgements.

J.E.P. gratefully acknowledges financial support from the Alexander-von-Humboldt Stiftung and the Spanish MEC (grant No. EX2001 11808094).


format: text

text:
# Probing Interstellar Dust Models Through SAXS (Small Angle X-Ray Scattering)

## 1 Introduction

The presence of dust in the general interstellar medium (ISM) is primarily inferred from the absorption, scattering and polarization of stellar and diffuse Galactic light, the infrared, submillimeter, and microwave emission spectrum from the diffuse ISM, the presence of broad solid state emission and absorption features in stellar and nebular spectra, and from the presence of X-ray scattering halos around X-ray sources. Indirect evidence for the presence of dust is inferred from the depletion of various refractory elements from the gas phase of the ISM, the presence of ”presolar” dust grains in the solar system, and the role of dust as a catalyst for chemical reactions, and as a provider of photoelectrons in the heating of the neutral ISM. Comprising less than one percent of the mass of the interstellar gas, dust nevertheless plays a major role in the radiative appearance of galaxies, and the thermal and chemical balance of their interstellar medium.

Ideally, an interstellar dust model should be derived by simultaneously fitting all observational constraints such as the average interstellar extinction, the diffuse IR emission, the UV to optical wavelength dependent albedo, polarization, and scattering phase function, and the intensity and spatial distribution of X-ray halos. It must consist of particles with known optical, physical, and chemical properties, and require no more than the ISM abundance of any given element to be locked up in the dust. In practice, most interstellar dust models have been constructed by deriving the abundances and size distributions of some well studied solids, such as graphite or silicates, from select observations such as the average interstellar extinction, the polarization, or the diffuse infrared emission, and checked for consistency with other observational constraints such as the wavelength dependent albedo and interstellar abundances.

The first widely used dust model derived in such manner consisted of spherical, bare silicate and graphite particles characterized by a $`\mathrm{d}n/\mathrm{d}aa^{3.5}`$ power law distribution in grain radius $`a`$ in, respectively, the 0.025 to 0.25 $`\mu `$m, and the 0.005 to 0.25 $`\mathrm{\mu m}`$ size intervals (Mathis, Rumple, & Nordsieck 1977; MRN). Using the dust optical constants of Draine & Lee (1984; DL84), the MRN model provides a good fit to the average interstellar extinction curve. Assuming solar composition for the ISM, the model requires essentially all the interstellar carbon, C/H = 370 ppm (parts per million), and all the magnesium, silicon, and iron, Mg/H, Si/H, Fe/H = {34, 35, 28}, to be locked up in dust.

The shortcomings of deriving an interstellar dust model from only one set of observational constraints were demonstrated by the results of the Infrared Astronomical Satellite (IRAS) all sky survey, which provided the average infrared (IR) emission spectrum at 12, 25, 60, and 100 $`\mu `$m from the diffuse ISM. The observations showed an excess of 12 and 25 $`\mu `$m emission over that expected from dust heated by the local interstellar radiation field (ISRF) and radiating at the equilibrium dust temperature. Draine & Anderson (1985) suggested that the MRN grain size distribution should be extended to very small grains (VSG) with radii of $``$ 5 Å which undergo temperature fluctuations when heated by the ISRF. Allamandola, Tielens, & Barker (1985), and Léger & Puget (1984) identified these VSG with polycyclic aromatic hydrocarbon (PAH) molecules whose presence in the ISM was inferred from the ubiquitous solid state emission features at 3.3, 6.7, 7.6, 8.6, and 11.3 $`\mu `$m. Other possible candidates for the carriers of these IR emission features were summarized in a review by Tokunaga (1997).

The first models that attempted to fit in a self–consistent manner the interstellar extinction as well as the diffuse IR emission using PAHs as an interstellar dust component were presented by Désert, Boulanger, & Puget (1990) and Siebenmorgen & Krügel (1992). However, their models used an ad hoc formula to characterize the optical and UV extinction of the PAHs.

Interstellar polarization provides additional constraints on interstellar dust particles (Kim & Martin 1995, 1996; Li & Greenberg 1997). In particular, the latter authors presented an interstellar dust model consisting of PAHs, and cylindrical silicates coated with an organic refractory mantle. Their model satisfies the interstellar extinction, polarization, and solar abundance constraints. It did not attempt to fit the diffuse IR emission, and used particles with hypothetical optical properties to represent the far-UV extinction and the 2175 Å extinction bump.

The Diffuse IR Background Experiment (DIRBE) and Far Infrared Absolute Spectrophotometer (FIRAS) instrument on board the Cosmic Background Explorer (COBE) satellite provided the most extensive wavelength coverage (3.5 to 1000 $`\mu `$m) of the IR emission from the diffuse ISM, and indirect evidence for the emission from PAHs from this phase of the medium (Dwek et al. 1997). Dwek et al. (1997) attempted to fit the interstellar extinction and diffuse IR emission using a mix of bare silicate and graphite particles with DL84 optical constants, and PAHs with the UV and optical properties of Désert, Boulanger, & Puget (1990). The model failed to reproduce the observed interstellar extinction, primarily due to the non-physical nature of the UV-optical properties adopted for the PAHs, and the adherence to a power law representation for the grain size distribution. Draine & Li (2001) and Li & Draine (2001; LD01) improved on this type of model by using a more realistic characterization of the optical properties of the PAHs, based on laboratory measurements. However, they used a predetermined functional form for the grain size distribution (Weingartner & Draine 2001; WD01), and did not include interstellar abundance constraints in the fitting procedure. Consequently, their model requires an excessive amount of Mg, Si, and Fe to be locked up in dust. The Si/H = Mg/H = Fe/H abundance in the dust is about 50 ppm (parts-per-million = 10-6), almost twice the amount of iron if the ISM has solar abundances.

Since interstellar abundances play a crucial role in constraining the amount of matter that can be locked up in the dust, Zubko, Dwek, & Arendt (2003; ZDA03) included interstellar abundances as an explicit constraint in the construction of dust models. Previously, models were just examined for consistency with abundance constraints, and non-compliances were given little weight because of the uncertainties in the determination of interstellar abundances. ZDA03 explored a variety of ISM abundances and potential dust compositions. The models are physical in the sense that they use measured optical constants or observed radiative properties to characterize the optical properties of the various dust constituents.

We first present, in §2, a brief summary of the ZDA03 dust models. An important result of their investigation is that there exist many, distinct, models that simultaneously satisfy the extinction, IR emission, and abundance constraints. X-ray halos provide an important additional constraint on vaible interstellar dust models. They are produced by the scattering of X-rays primarily off large dust grains, a region of grain sizes largely inaccessible at other wavelengths. In §3 we briefly summarize the basic X-ray scattering properties of individual dust grains, and the collective effects of many grains along the line of sight to an X-ray source on the intensity and spatial profile of its scattering halo. We review recent modeling efforts of the halo around Nova Cygni 1992, and then present preliminary models for the halo around the X-ray source GX13+1, showing that is a powerful discriminator between different classes of dust models, significantly narrowing the number of viable models presented by ZDA03 (§4). A brief summary of the paper is presented in §5.

## 2 Interstellar Dust Models Consistent With Extinction, IR Emission, and Abundance Constraints

Given the composition and physical and optical properties of the dust particles, the problem of simultaneously fitting the model to a set of observational constraints is a typical ill-posed inversion problem, in which the grain size distribution is the unknown. ZDA03 solved the inversion problem using the method of regularization, and considered five different dust compositions as potential model ingredients: (1) PAHs; (2) graphite; (3) hydrogenated amorphous carbon of type ACH2; (4) silicates (MgSiFeO4); and (5) composite particles containing different proportions of silicates, organic refractory material (C8H8O4N), water ice (H2O), and voids.

These five dust compositions were combined to create five different classes of dust models.
The first class consists of PAHs, and bare graphite and silicate grains, and is identical to the carbonaceous/silicate model recently proposed by LD01.
The second class of models contains composite particles in addition to PAHs, bare graphite and silicate grains.
The third and fourth classes of models comprise the first and second classes but with the graphite grains completely replaced by amorphous carbon grains.
In the fifth class of models the only carbon is in PAHs and in the organic refractory material in composite grains. That is, the model comprises only PAHs, bare silicate, and composite particles.

Models including only PAHs and bare grains were designated as BARE, and those containing additional composite particles as COMP. The BARE and COMP models were further subdivided according to the composition of the bare carbon particles used in the model: graphite, amorphous carbon, and no carbon, designated respectively by -GR, -AC, and -NC. For all model classes, ZDA03 performed the fit using three different sets of ISM abundances: solar, B stars, and F and G stars (correspondingly designated as -S, -B, and -FG submodels). So, for example, a BARE model with graphite grains derived by assuming the solar ISM abundances is designated as BARE-GR-S, or a COMP model with no carbon obtained by assuming B star ISM abundances is designated as COMP-NC-B.

The method and the results of the analysis are shown in the form of a work flow chart in the top half of Figure 1 which is an edited reproduction of the poster exhibited at this conference.

The left part of the Work Flow Chart depicts the observational constraints used in the work of ZDA03. They include the average interstellar extinction, the diffuse IR emission, and the abundance constraints.

The right part of the Work Flow Chart summarizes the different dust compositions used in the model.

The center of the Work Flow Chart depicts the method of regularization, and the different classes of dust models used in the fit to the observational constraints.

Model results are depicted below the Work Flow Chart. Of particular interest are the curves in the upper row, showing the derived grain size distributions for the different models. The lowest row shows the abundances of the different elements locked up in the dust. B stars abundances place the most stringent constraint on the amount of carbon that can be locked up in the dust, less than $``$ 160 ppm. A model that could fit this constraint (as well as the extinction and IR emission) consisted of PAH, bare silicate, and composite grains (no bare carbon particles!). In this model (model COMP-NC-B) the 2175 Å extinction bump was entirely produced by PAHs.

Figure 2 offers a detailed comparison of the PAH and bare graphite and silicate model derived by LD01 and ZDA03. The main peak in the LD01 graphite size distribution peaks at larger grain sizes ($`0.3\mu `$m), whereas the ZDA03 peaks at $`0.1\mu `$m. The ZDA03 size distribution is smoother, only double–peaked: one peak corresponding to graphite grains, the other to PAHs. It lacks the more complex triple–peak distribution of the LD01 model. Both models consume about equal amounts of carbon in PAHs and graphite with C/H $``$ 240 ppm. The models differ significantly in the size distribution of the silicate particles. The LD01 silicate size distribution is steeper and has a larger amplitude, whereas the ZDA03 distribution has a slope similar to the MRN size distribution, a significantly lower peak, and a larger mass fraction of small silicate grains. The amount of Si/H, Mg/H, and Fe/H in the ZDA03 model is about 33 ppm for each element, consistent with solar abundances, whereas the LD01 model requires an abundance of $``$ 50 ppm for each element, almost twice the available iron abundance if the ISM had solar composition.

The results summarized above show that there are many classes of interstellar dust models that provide good simultaneous fits to the far-UV to near-IR extinction, thermal IR emission, and elemental abundance constraints. The models can be grouped into two major categories: BARE and COMP models. The latter are distinguished from the former by the fact that they contain a population of composite particles which generally comprise of larger grains than bare particles. Furthermore, a significant fraction of the volume of composite grains consists of voids, so their effective electron density is significantly smaller than that of the bare particles. Since X-ray halos are primarily produced by scattering off electrons in large grains, they may provide an excellent discriminator between viable dust models.

## 3 The Production of X-Ray Halos

The use of X-ray halos as diagnostic tool for determining dust properties, the morphology of the ISM, and distances to X-ray sources, have been previously discussed in the literature by, among others, Overbeck (1965), Stecher & Williams (1966), Trümper & Schönfelder (1973), Alcock & Hatchett (1978), Mauche & Gorenstein (1986), Mathis & Lee (1991), Predehl & Schmitt (1995), Predehl & Klose (1996), and Predehl et al. (2000). Here we will concentrate on the two main factors that determine the X-ray halo profile and intensity: (1) the scattering properties of the individual dust grains, and (2) the collective effects of the many individual dust grains along the line of sight (LOS) to the X-ray source.

### 3.1 The Scattering of X-Rays by Single Dust Particles

The scattering of X-rays is usually described by the Rayleigh-Gans theory which is valid under the following two assumptions (van der Hulst 1981):
(1) $`|m1|1`$, where $`m`$ is the refractive index of the dust particle. This conditions ensures that the incident waves are not significantly reflected off the dust; and (2) $`\mathrm{\Delta }\mathrm{\Phi }=(4\pi a/\lambda )|m1|1`$, ensuring that there is no appreciable phase shift across the grain, so that each element in the grain is exposed to the same phase and amplitude of the incident wave.

The second condition is the more stringent one which can be written as (Alcock & Hatchett 1978):

$$E_x(keV)3a(\mathrm{\mu m})\left(\frac{\rho }{3\mathrm{g}\mathrm{cm}^3}\right),$$
(1)
where $`E_x`$ is energy of the incident X-rays. The Rayleigh–Gans approximation is therefore valid for X-ray energies in excess of a few keV, and dust particles with radii $`a<1\mathrm{\mu m}`$.

When these conditions are met, each electron in the grain can be considered as an independent scatterer. Waves scattered in a given direction interfere with each other because of the different positions of the scatterers within the grain. The scattering is strongly peaked in the forward direction since, when $`\theta \lambda /a`$, the phase differences between the scattering electrons across the grain are very small.

The interference effects can be readily calculated for simple grain geometries. Consider an X-ray incident on a dust particle, which for illustrative purposes is taken to be a sphere of radius $`a`$ (see Figure 3). Since X-ray energies are in general much larger than the electron binding energy to an atom, the electrons can be considered as free particles, confined to the volume of the dust grain, and vibrating at the frequency of the incident X-ray. Let $`N`$ be the position in the dust of an electron scattering the X-ray at an angle $`\theta `$ with respect to the direction of the incident beam. The path difference between the scattered rays eminating from positions $`N`$ and $`P`$ is given by $`\overline{OP}\overline{MN}=\underset{¯}{r}(\underset{¯}{k}\underset{¯}{k}_0)\underset{¯}{K}\underset{¯}{r}`$, where $`\underset{¯}{k}_0`$ and $`\underset{¯}{k}`$ are unit vectors in the direction of the incident and scattered beam, respectively. From Fig. 3 we see that $`|\underset{¯}{K}|=2\mathrm{sin}(\theta /2)`$.
The phase difference between the two rays is therefore given by

$`\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{2\pi }{\lambda }}\underset{¯}{K}\underset{¯}{r}`$
$`=`$ $`{\displaystyle \frac{4\pi r}{\lambda }}\mathrm{sin}\left({\displaystyle \frac{\theta }{2}}\right)\mathrm{cos}\alpha `$
where $`\alpha `$ is the angle between $`\underset{¯}{r}`$ and $`\underset{¯}{K}`$. The scattering is coherent, so the amplitudes are first added and then squared to get the intensity of the scattered light. The scattering factor, defined as the normalized sum of the scattering amplitudes in a given direction, is given by:

$`f(\theta )`$ $`=`$ $`{\displaystyle \rho (r)\mathrm{exp}(i\mathrm{\Phi })dV}`$
$`=`$ $`{\displaystyle \rho (r)\mathrm{exp}\left[i\kappa r\mathrm{cos}\alpha \right]2\pi r^2dr\mathrm{sin}\alpha \mathrm{d}\alpha }`$
where $`\kappa 4\pi \mathrm{sin}(\theta /2)/\lambda `$, d$`V=2\pi r^2`$d$`r\mathrm{sin}\alpha `$d$`\alpha `$ is the differential volume element in a spherical coordinate system in which the $`z`$axis is in the $`\underset{¯}{K}`$ direction, and $`\rho (r)`$ is the electron density in the solid. For a homogeneous dust particle the scattering factor simplifies to (see van de Hulst 1981, §7.21):

$`f(\theta )`$ $`=`$ $`{\displaystyle 4\pi \rho (r)r^2\frac{\mathrm{sin}(\kappa r)}{\kappa r}dr}`$
$`=`$ $`{\displaystyle \frac{3N_e}{x^3}}\left[\mathrm{sin}xx\mathrm{cos}x\right]=3N_e\left({\displaystyle \frac{\pi }{2x^3}}\right)^{1/2}J_{3/2}(x)`$
where $`x=\kappa a`$, and $`J_{3/2}(x)`$ is the Bessel function of order $`n=3/2`$. The scattering amplitude is proportional to $`f^2(\theta )N_e^2`$, where $`N_e`$ the total number of free electrons in the grain. The scattering amplitude for the sphere is the 3-dimensional analog of the scattering amplitude for a circular aperture (Airy disk) which has a diffraction pattern proportional to $`J_1^2(x)/x^2.`$

When the assumptions of the Rayleigh-Gans approximation are not fulfilled, for example for X-rays emitted in the ROSAT $`\frac{1}{4}`$ and $`\frac{3}{4}`$ keV energy band, the scattering must be calculated using detailed Mie calculations, which are valid even for particles with $`a\lambda `$ (Smith & Dwek 1998). The scattering cross section exhibits then a less than quadratic dependence on $`N_e`$.

Regardless of the scattering regime and the numerical procedure used in the calculations, the scattering amplitude is always strongly peaked in the forward direction. This effect is shown in Figure 4, which depicts the polar scattering diagram for 2.9 keV X-rays incident on 0.001, 0.01, 0.1, and 1.0 $`\mathrm{\mu m}`$ silicate grains. Figure 5 shows the dependency of the scattering phase function, defined as $`\sigma (\theta )/\pi a^2`$, versus scattering angle for 0.01 and 0.1 $`\mathrm{\mu m}`$ graphite grains. The figure shows that the scattering phase function is brightest at $`\theta `$ = 0, steeply declining at angles $`>\lambda /a`$.

### 3.2 The Collective Effects of Dust Particles on the Halo Profile and Intensity

Having summarized the size dependent scattering properties of individual dust particles, we will now examine the collective effect of a population of dust particles along the LOS to the source. The intensity and spatial profile of X-ray halo will depend microscopically on the size distribution and properties of the individual dust grains, and macroscopically on the morphology of the ISM.

#### 3.2.1 Clumpy ISM.

Consider first the formation of a halo in a clumpy ISM. Figure 6a compares the scattering geometry for two single dust clouds located at distances $`d_1`$ and $`d_2`$ in front of an X-ray source $`S`$ as viewed by an observer located at distance $`D`$ from the source. For sake of clarity, the scattering angles $`\theta _1`$ and $`\theta _2`$ are highly exaggerated. The X-ray flux incident on Cloud 1 is larger than that incident on Cloud 2 ($`d_1<d_2`$), so Cloud 1 will have a higher surface brightness than Cloud 2. However, the scattering angle $`\theta _1`$ is larger than $`\theta _2`$, which can further affect the relative surface brightness of the clouds. To illustrate this effect we will assume that the population of dust particles in each of the two clouds has a bimodal size distribution sharply peaked at radii $`a`$ = 0.01 and 0.1 $`\mathrm{\mu m}`$. If Cloud 1 is sufficiently close to the LOS so that $`\theta _1<200\mathrm{}`$, then both size populations in each cloud will contribute to the scattering (see Figure 5). However, if $`\theta _1>200\mathrm{}`$, the contribution of the population of 0.1 $`\mathrm{\mu m}`$ dust particles to the scattering from Cloud 1 will drop dramatically. However, the same population of dust particles may still contribute to the scattering intensity from Cloud 2, provided that $`\theta _2<200\mathrm{}`$. This simplified scenario illustrates the fact that: in a clumpy ISM, one can control the intensity and spatial profile of a halo by adjusting the grain size distribution and the distance of the dust from the source.

#### 3.2.2 Uniform ISM.

In an ISM in which the dust is uniformly distributed along the LOS, all scattering angles are possible (see Figure 6b). Without the flexibility of biasing the scattering to a set of discrete angles, one can now only affect the intensity of the halo at, for example, small scattering angles by changing the abundance of small angle scatterers, that is, by adjusting the population of large dust particles.

#### 3.2.3 Summary of Factors Determining the Halo Intensity and Profile.

To summarize, the main factors that play a role in the formation of a halo are:
(1) Source properties: The energy spectrum of the X-ray source affects the halo properties in broad–band X-ray observations. Source variability will cause correlated spatial variations in the intensity and profile, allowing for distance determination and dust distribution to the X-ray source (Alcock & Hatchett 1978, Predehl et al. 2000).
(2) ISM properties: The distribution of the dust along the LOS, and the total dust column density constrain the scattering geometry.
(3) Dust properties: The dust composition (heavy elements can bind electrons more strongly that light ones - a weak effect), the density of the dust grains, and the grain size distribution are the most fundamental factors determining the halo profile.

## 4 Constraining Dust Models With X-Ray Halo Observations

### 4.1 Nova Cygni 1992

Nova Cygni 1992 was observed with the ROSAT satellite, which revealed the existence of an extensive halo surrounding a bright source of soft X-rays in the 0.4 keV energy band. The interstellar reddening to this source has been estimated to be between $`E(BV)`$ = 0.19 and 0.31, corresponding to an H-column density of $`N_\mathrm{H}`$ = 1.1 and 1.8 $`\times 10^{21}`$ cm-2, respectively (Mathis et al. 1995). The source coordinates are {l, b} = {89.13$`\mathrm{deg}`$, 7.82$`\mathrm{deg}`$}, and its most accepted distance is about 3.2 kpc. The reddening provides an integral constraint on the grain size distribution weighted towards small grains, whereas the X-ray scattering can constrain the large end of the grain size distribution, where $`E(BV)`$ is largely flat. Simultaneously fitting the observed reddening, halo intensity, and abundance constraints is a delicate balancing act of having enough dust to produce the reddening without overproducing the halo.

Mathis et al. (1995) found that reproducing the observed optical extinction without overproducing the halo intensity, requires the dust to be composed of silicate and amorphous carbon, with voids making up a significant fraction of the volume of the large ($`>0.1\mathrm{\mu m}`$) grains. Fluffy grains scatter X-rays less efficiently than solid grains of the same size, and can still produce the observed optical extinction.

Smith & Dwek (1998) challenged the result that fluffy grains are required to simultaneously fit the extinction and halo towards Nova Cygni 1992. Mathis et al. (1995) used the Rayleigh-Gans approximation to calculate the X-ray scattering. This approximation breaks down at the energy of the ROSAT observations. Detailed Mie calculations show that the X-ray scattering is less efficient than that predicted by the Rayleigh-Gans approximation, and Smith & Dwek succeeded in reproducing the reddening and the general halo profile with an MRN mixture of bare silicate and graphite grains.

Observational evidence for the presence of large ($`>1\mathrm{\mu m}`$) interstellar dust grains entering the solar system, summarized by Frisch et al. (1999), prompted Witt, Smith, & Dwek (2001) to reanalyze the Nova Cygni 1992 halo in search for evidence for the prevalence of such large dust particles in the general ISM. Adopting a uniform distribution of dust along the LOS, an H-column density of $`2.1\times 10^{21}`$ cm-2, and a reddening value of $`E(BV)`$ = 0.36, they found that an MRN size distribution, extended to radii larger than 2 $`\mathrm{\mu m}`$ could fit the X-ray halo without violating solar abundance constraints. Their fit was more detailed than the one by Smith & Dwek (1998), and driven by the halo intensity at the smallest scattering angle in their data of $`100\mathrm{}`$. Most of the dust mass in the extended MRN distribution, XMRN, is concentrated in the large grains, enhancing the halo intensity at small scattering angles while weakening it at the larger ones. The problem with the XMRN distribution is that it produced a large visual optical depth, leading to an extinction curve with $`R_VA_V/E(BV)=6.1`$, twice the value for the average ISM extinction of 3.1. This would suggest a peculiar size distribution along the Nova Cygni 1992 LOS, different from the average ISM. Equivalently, the large extinction produced an $`A_V/N_\mathrm{H}`$ value of $`1\times 10^{22}`$ cm2, twice the canonical ISM value.

Draine & Tan (2003) reanalyzed the X-ray halo around Nova Cygni 1992 for six different epochs. They presented halo data down to a scattering angle of $`55\mathrm{}`$, providing stronger constraints on the population of large dust particles, or the LOS distribution of the dust particles. Standard LD01 dust uniformly distributed along the LOS underproduced the halo at small scattering angles, while overproducing it at larger ones (see Figure 15 in their paper). The halo intensity could be increased at small scattering angles by the addition of larger grains in the WD01 size distribution without affecting the halo intensity at large scattering angles. However, they preferred the dust model consisting of LD03 dust with an WD01 size distribution in which 30% of the dust is in one cloud located close to the source at 0.05 of the source distance to the observer. The model could reproduce the average Galactic extinction with $`R_V`$ = 3.1, provided that $`E(BV)0.20`$, leading to a H-column density of $`N_\mathrm{H}1.1\times 10^{21}`$ cm-2 for a canonical $`A_V/N_\mathrm{H}`$ value of $`5.5\times 10^{22}`$ cm2. However, such a cloud would have to be at a height of $``$ 410 pc above the plane for a source distance of 3.2 kpc. This height is significantly larger than the $``$ 130 pc scale height of the Galactic dust layer (Cox 1999). Draine & Tan (2003) place the source at a shorter distance of 2.1 kpc, which will still put the cloud at a height of $``$ 270 pc above the plane.

### 4.2 The X-Ray Binary GX 13+1

The high angular resolution of the Chandra X-ray telescope makes it an ideal tool for studying X-ray halos. Smith, Edgar, & Shafer (2002, SES02) presented halo observations of GX 13+1 at energies of 2.1, 2.9, and 3.7 keV. In addition to the halo properties, an important constraint on dust models is the total H-column density to the source, estimated to be $`N_\mathrm{H}=2.9\pm 0.1\times 10^{22}`$ cm-2 (see SES02 for references). SES02 found that the halo could be fitted by an MRN or a WD01 grain size distribution. The observations ruled out fluffy dust models, as well as the XMRN model used by Witt, Smith, & Dwek (2001).

The Chandra observations of GX 13+1 provide a useful means for discriminating between the different dust models of ZDA03. Figure 7 shows the fit to the 2.9 keV halo intensity profile for a dust model consisting of PAHs and bare silicate and graphite grains (BARE-GR-S), uniformly distributed along the LOS. Similarly good fits were obtained with the same dust model for the halo intensity profile at 2.1 and 3.7 keV.

Figure 7 also shows that the composite dust model COMP-AC-S provides an equally good fit as the BARE-GR-S model to the halo profile. However, it requires a H-column density of $`4.1\times 10^{22}`$ cm-2, significantly above the observational value. The reason for this large column density can be understood by examining the scattering phase function of the different dust constituents in each of the two dust model (Figure 8). COMP dust models contain a population of composite particles which generally comprise larger grains than bare particles. Their scattering phase function exhibits a prominent rise at small scattering angles so that the halo intensity is dominated by composite particles for $`\theta <100\mathrm{}`$. However, a significant fraction of their volume contains voids, so that their scattering efficiency is significantly reduced compared to same–sized bare particles. In order to fit the halo intensity with such a population of inefficient scatterers their total abundance needs to be increased. However, the ISM abundance constraints prohibit such change, so the only way of fitting the halo with composite particles without violating the abundance constraints is to increase the H-column density as well. But this quantity is also constrained, by X-ray observations. The overall result is that the halo around GX 13+1 cannot be reproduced by COMP models. This example illustrates the need for reliable H-column density measurement in order to discriminate between otherwise viable dust models. Color excess ratios toward the X-ray sources will also be extremely useful additional constraints.

## 5 Summary

Dust models depend on the observational constraints that are explicitly included in the fitting procedure, and whether the procedure constrains the grain size distribution to be of a given functional form. Recently, Zubko, Dwek, & Arendt (2003) have constructed dust models explicitly constrained by the average interstellar extinction, diffuse IR emission, and ISM abundance constraints, without an a priori functional form for the grain size distribution. Their results show that these constraints are insufficient for identifying a unique dust model. The radial intensity distribution and spectrum of X-ray scattering halos provide additional constraints that help characterize the composition and relative numbers of the larger dust grains. However, since the extinction, IR emission, and abundance constraints represent averages of many different lines of sight, a large observational data base of X-ray halos is needed to render X-ray scattering an equally relevant constraint.

Acknowledgments
We thank Adolf Witt and the Scientific Organizing Committee for creating a very stimulating conference in a very beautiful and invigorating setting. This work was supported by NASA’s Astrophysics Theory Program NRA 99-OSS-01.


text:
# 1 Introduction

## 1 Introduction

A gas of ground state 3He atoms in which a high degree of nuclear polarization is achieved offers an incredibly rich playground in various fields of science, from statistical or nuclear physics to biophysics and medicine . Depending on the targeted application, the degree of nuclear polarization, the sample density, or the production rate of polarized atoms should be optimized. A recent application, which may have an important impact on the diagnosis of pulmonary diseases, is polarized gas magnetic resonance imaging (MRI) . Clinical studies to demonstrate the relevance of this new tool are under way in Europe and in the United States. Yet, if a wide expertise exists in MRI to adapt the existing imaging techniques to the case of polarized gases, the gas preparation remains a critical stage to be transferred from physics laboratories to hospitals. Two methods are presently used to polarize 3He: spin-exchange with optically pumped alkali atoms  and pure-He metastability exchange optical pumping (MEOP) . In standard conditions, MEOP is performed at low pressures ($`1`$ mbar) in a guiding magnetic field of the order of $`1`$ mT. Circularly polarized light at $`1083`$ nm, corresponding to the $`2^3`$S-$`2^3`$P transition of 3He, is used to transfer angular momentum to the atoms and nuclear polarization is created by hyperfine coupling in the metastable $`2^3`$S state. Through metastability exchange collisions, nuclear polarization builds up in the ground state. The steady-state nuclear polarization obtained by MEOP in standard conditions rapidly decreases if the pressure of the sample exceeds a few mbar (see below, Fig.5-a) . Therefore a delicate polarization-preserving compression stage is necessary for MRI where the gas should be at atmospheric pressure for inhalation, and for all applications needing a dense sample. In this letter, the MEOP scheme is shown to withstand large hyperfine decoupling. A strong magnetic field of $`1.5`$ T actually improves its performances with respect to standard low-field optical pumping. At $`1.33`$ mbar, high nuclear polarizations of the order of $`80\%`$ are routinely obtained with much lighter experimental constraints. At higher pressures, the achieved nuclear polarizations are dramatically increased compared to published low-field results. An elementary model with simple rate equations is used to account for these results.

## 2 Experimental setup and methods

Experiments are performed in the bore of the $`1.5`$ T superconducting magnet of a clinical MRI system. The experimental apparatus is sketched in Fig.1.

The helium gas is enclosed in a sealed cylindrical Pyrex cell (diameter $`5`$ cm, length $`5`$ cm). Four cells filled with $`1.33`$ mbar, $`8`$ mbar, $`32`$ mbar and $`67`$ mbar of pure 3He are used. A radio-frequency discharge at $`3`$ MHz is sustained in the gas by external electrodes, leading to metastable atom densities $`n_m`$ in the $`0.3`$-$`3\times 10^{10}`$atoms/cm3 range, depending on the applied voltage and on the gas pressure. The optical pumping laser is either a $`50`$ mW single mode laser diode amplified to $`0.5`$ W , or a broadband fiber laser ($`1.63`$ GHz FWHM) delivering $`2`$ W . The pump beam is back-reflected to enhance its absorption, which is monitored on the transmitted beam with a photodiode. At the entrance of the cell, the Gaussian transverse intensity profile of the pump beam has a FWHM of the order of $`2`$ cm. A weak probe beam from a single mode laser diode is used to measure the nuclear polarization. It is linearly polarized perpendicularly to the magnetic field ($`\sigma `$ polarization). The discharge intensity is modulated at $`133`$ Hz, and the probe absorption is measured with a lock-in amplifier. Laser sources and electronics remain several meters away from the magnet bore, in a low-field region.

At $`1.5`$ T, due to Zeeman splitting, the energy levels of the $`2^3`$S and $`2^3`$P states are spread over $`80`$ and $`160`$ GHz respectively (Fig.2-a).

Hyperfine decoupling in the $`2^3`$S state is significant, so that the eigenstates of the Hamiltonian show only little mixing between different eigenstates $`|m_J,m_I`$ of the decoupled spin system, where $`m_J`$, $`m_I`$, and $`m_F`$ denote the magnetic quantum numbers for the electronic, nuclear, and total angular momentum, respectively. As shown in Fig.2-a, the $`2^3`$S sublevels are arranged into three pairs of quasi-degenerate levels of increasing energies (A1,A2), (A3,A4), and (A5,A6) that correspond respectively to $`m_J`$=$`1`$, $`0`$, and $`1`$ in the completely decoupled limit $`B\mathrm{}`$. For more details about the $`2^3`$S level structure and the analytical expressions of eigenstates and energies, we refer the reader to . The absorption spectra at low magnetic field and at $`1.5`$ T are displayed in Figs.2-b and c, respectively. In standard MEOP, very high nuclear polarizations are obtained using C8 or C9 lines . Comparable polarizations are achieved at $`1.5`$ T using the $`\sigma ^{}`$-strong pump line displayed in Fig.2-c. All the results presented in this work are obtained with this pump transition. The performances and efficiencies of other optical pumping transitions at $`1.5`$ T will be reported elsewhere. The pump simultaneously addresses the four $`2^3`$S sublevels A1 to A4. Population transfer into the pair (A5,A6) is achieved by the following sequence : laser excitation, collisional redistribution in the $`2^3`$P state and spontaneous emission. The ground state nuclear polarization $`M`$ is defined as $`M`$=$`(n_+n_{})/(n_++n_{})`$ where $`n_+`$ and $`n_{}`$ denote populations of the $`m_I`$=$`1/2`$ and $`m_I`$=$`1/2`$ nuclear spin states, respectively. In the absence of optical pumping, metastability exchange collisions impose a spin temperature distribution for the $`2^3`$S sublevel populations, proportionally to $`e^{\beta m_F}`$ where $`e^\beta `$=$`n_+`$/$`n_{}`$=$`(1+M)/(1M)`$ .

The populations of sublevels A5 and A6, not addressed by the pump, are probed to measure $`M`$. Examples of probe absorption spectra for an unpolarized and an optically-pumped steady-state situation are shown in Fig.3-a. $`M`$ is inferred from the relative heights of the absorption peaks. The build-up of the polarization, as well as its decay when the pump is turned off, are monitored by tuning the probe laser frequency on the probe transition starting from the A5 ($`m_F`$=$`3/2`$) sublevel (Fig.3-b). These measurement procedures operate at arbitrary magnetic field and pressure .

## 3 Results

The steady-state polarization $`M_{eq}`$ and the polarization build-up time constant $`t_b`$ in the $`1.33`$ mbar cell are shown in Figs.4-a and b as a function of the discharge-induced decay time $`T_1`$.

Over a wide range of moderate to weak discharges ($`T_1`$ ranging from $`300`$ s to $`1500`$ s), $`t_b`$ (ranging from $`60`$ to $`350`$ s) is proportional to $`T_1`$ and the polarization achieved with the broadband $`2`$ W laser is high, about $`80\%`$, independently of $`T_1`$. This behavior is specific to the high-field optical pumping, and contrasts with the standard low-field situation where a very weak discharge is required to obtain such large nuclear polarizations. For the strongest discharges, build-up times decrease ($`t_b`$ ranging from $`15`$ to $`55`$ s) and steady-state polarizations are lower. Fig.4-c and d show the influence of the pump laser power for a given discharge ($`T_1`$=$`300`$ s). A laser power as low as $`0.5`$ W is sufficient for the polarization and the build-up time to almost reach their asymptotic values.

Similar studies of MEOP have been performed at higher pressures, where $`T_1`$ between $`300`$ and $`2600`$ s are measured. Selected results for a weak and a strong discharge at full laser power are shown in Fig.5 together with published low-field results. The polarizations obtained at high pressures are, to our knowledge, record MEOP values. The polarization build-up times weakly depend on 3He pressure, in contrast with low-field MEOP .

## 4 Discussion

An extension of the detailed model for standard MEOP  to high-field conditions  is required to compute the populations of all atomic sublevels. Here, for simplicity, an elementary model is used to account for the main observed features. We assume that (i) atoms are fully pumped into the (A5,A6) pair, and (ii) the populations of sublevels not addressed by the pump laser are imposed by the ground state spin temperature which only depends on $`M`$: $`a_5`$=$`(1+M)/2`$ and $`a_6`$=$`(1M)/2`$. The sublevel A5 is totally oriented ($`m_J`$=1, $`m_I`$=1/2) and carries a nuclear angular momentum $`I_z`$=$`\mathrm{}/2`$, while A6 has a small component on ($`m_J`$=0, $`m_I`$=1/2) and a large component on ($`m_J`$=1, $`m_I`$=-1/2) and thus carries a nuclear angular momentum $`I_z`$=$`\mathrm{}(ϵ1)/2`$ with $`ϵ`$=$`1\times 10^2`$ at $`1.5`$ T . The rate equation for $`M`$, resulting from relaxation and metastability exchange, is then:

$$\frac{dM}{dt}=\frac{2I_z/\mathrm{}M}{T_e}\frac{M}{T_1}\text{with}I_z=\frac{\mathrm{}}{2}(M+ϵ\frac{1M}{2}),$$

where $`1/T_e`$ is the metastability exchange collision rate for a 3He atom in the ground state ($`1/T_e`$=$`n_m\times 1.53\times 10^{10}`$cm3/s), and $`2I_z/\mathrm{}`$ is the nuclear polarization in the $`2^3`$S state. One infers a steady-state polarization $`M_{eq}`$=$`(1+2T_e/(ϵT_1))^1`$ and a build-up time $`t_b`$=$`2T_eM_{eq}/ϵ`$. Using values of $`n_m`$ and $`T_1`$ measured in the absence of pumping beam, the predicted polarization is $`M_{eq}`$$``$$`1`$, at all pressures. The build-up times are in the range $`20`$-$`300`$ s for the low pressure cell, and in the range $`15`$-$`40`$ s for the three high pressure cells. Although this elementary model is clearly not sufficient to predict $`M_{eq}`$, it accounts reasonably well for the observed dynamics. Its domain of validity and accuracy are estimated from detailed rate equations for the six $`2^3`$S-sublevel populations. We find that in our experimental conditions and for the observed range of nuclear polarization, $`2I_z/\mathrm{}`$ given by the simple model differs from the exact value by a factor not exceeding 2, depending on $`M`$ and on the gas pressure. This difference results from incomplete population transfer into (A5,A6) as well as from deviations of the order of $`ϵ`$ of the populations $`a_5`$ and $`a_6`$ from their assumed spin-temperature values. In spite of its simplicity, the model sheds light on two key features: (i) The dynamics of optical pumping at $`1.5`$ T is intrinsically limited by hyperfine decoupling. (ii) The build-up time, at least in the explored range of parameters, weakly depends on pressure and is affected only through changes of $`n_m`$ and $`T_1`$.

For application purposes, production rates of polarized atoms per unit volume $`R_a`$=$`PM_{eq}/t_b`$ are compared to published results for standard MEOP conditions and similar sealed cells in Table 1.

At low pressure, production rates at high field are lower than those obtained with low-field optical pumping. Nevertheless, one can take advantage of the weak pressure dependence of $`M_{eq}`$ and $`t_b`$ at $`1.5`$ T to efficiently perform MEOP at higher pressure. By increasing the pressure from $`1.33`$ to $`32`$ mbar, a factor of $`10`$ in $`R_a`$ is gained and good production rates are recovered. For instance, gas in a $`250`$ cc cell at $`32`$ mbar can be polarized at $`40\%`$ within $`3`$ minutes. This amount of gas is suitable for small animal lung imaging after compression to atmospheric pressure. For human lung MRI, considerable scaling-up or accumulation of polarized gas remains necessary. However, optical pumping around $`50`$ mbar would considerably simplify the compression stage by reducing the compression ratio from 1:1000 down to 1:20.

An intrinsic advantage of the high-field MEOP scheme is that, due to the large Zeeman splittings in the $`2^3`$S-$`2^3`$P transition, the magnetic sublevels involved in optical pumping are selected by the frequency of the light, and not only by its polarization. High-field MEOP is therefore extremely robust against polarization impurity of the pumping light. This is a crucial issue for massive production of polarized 3He using high laser power, since imperfect light polarization can severely limit achieved polarizations at low field .

## 5 Perspectives

The nuclear polarization improvement observed at $`1.5`$ T for high pressures is plausibly due to the inhibition by hyperfine decoupling of relaxation channels in atomic and/or molecular states in the plasma, as suggested by preliminary results at $`0.1`$ T . Further experiments at different magnetic field intensities are planned to confirm this hypothesis. In this perspective, the present study provides a first set of data showing that, in spite of the large hyperfine decoupling in the $`2^3`$S state, MEOP at high field (i) still yields high nuclear polarizations at low pressures and (ii) extends the domain of its applicability to higher pressures, providing fair polarizations and high production rates. From a practical point of view, and in the perspective of a large scale medical use of polarized gases, the development of a 3He polarizer operating at $`1.5`$ T (a widely used magnetic field in MRI), and at tens of mbar (for simplified compression) could be an attractive choice.


text:
# Chapter 1 Cosmic ray acceleration at parallel relativistic shocks in the presence of finite-amplitude magnetic field perturbations

## Chapter 1 Cosmic ray acceleration at parallel relativistic shocks in the presence of finite-amplitude magnetic field perturbations

### Abstract

We study the first-order Fermi acceleration process at parallel shock waves by means of Monte Carlo simulations. A ’realistic’ model of the magnetic field turbulence is applied involving sinusoidal waves imposed on the background mean field parallel to the shock normal. The finite-amplitude magnetic field perturbations lead to locally oblique field configurations at the shock and the respective magnetic field compression. It results in modification of the particle acceleration process introducing some features observed in oblique shocks, e.g. particle reflections from the shock. For the parallel mildly relativistic shocks we demonstrate for the first time a (non-monotonic) variation of the accelerated particle spectral index with the turbulence amplitude.

### 1. Introduction

As discussed by Ostrowski  for nonrelativistic shocks, the presence of finite-amplitude magnetic field perturbations modify character of the diffusive particle acceleration at the shock wave with the mean field parallel to the shock normal. The effect arises due to locally oblique field configurations formed by long-wave perturbations at the shock front and the respective magnetic field compressions. As a result the mean particle energy gains may increase and the particles reflected from the shock front may occur. The same phenomena should occur at relativistic shocks \[cf. 10\].

In the simplified numerical simulations of the first-order Fermi acceleration at parallel mildly relativistic shocks the acceleration time scale reduces with increasing turbulence level, but no spectral index variation occurs \[1, see 2 for ultrarelativistic shocks\]. However, the considered acceleration models apply very simple modeling of the perturbed magnetic field effects by introducing particle pitch-angle scattering. The purpose of the present work is to simulate the first-order Fermi acceleration process at mildly relativistic shock waves propagating in more realistic perturbed magnetic fields, including a wide wave vector range turbulence with the power-law spectrum. The magnetic field is continuous across the shock, according to the respective jump conditions. This feature leads to substantial modifications of the acceleration process at parallel shocks.

Below the upstream (downstream) quantities are labeled with the index ‘1’ (‘2’).

### 2. Simulations

In the simulations trajectories of ultrarelativistic test particles are derived by integrating their equations of motion in the perturbed magnetic field \[for details see: Niemiec, Ostrowski, in preparation\]. A relativistic shock wave is modeled as a plane discontinuity propagating in electron-proton plasma. The magnetic field is defined upstream of the shock. It consists of the uniform component, $`\stackrel{}{B}_{0,1}`$, parallel to the shock normal and finite-amplitude perturbations imposed upon it. The perturbations are modeled as a superposition of 294 sinusoidal static waves of finite amplitudes \[cf. 10\]. They have either a flat $`(F(k)k^1)`$ or a Kolmogorov $`(F(k)k^{5/3})`$ wave power spectrum in the (wide) wave vector range $`(k_{min},k_{max})`$, where $`k_{max}/k_{min}=10^5`$. The shock moves with the velocity $`u_1`$ with respect to the upstream plasma. The downstream flow velocity $`u_2`$ and the magnetic field structure are obtained from the hydrodynamic shock jump conditions. Derivation of the shock compression ratio as measured in the shock rest frame, $`R=u_1/u_2`$, is based on the approximate formulae derived in Ref. . In the analysis of the acceleration process the particle radiative (or other) losses are neglected.

### 3. Results for parallel shocks

In Fig. 1 we present particle spectra for the parallel shock wave with $`u_1=0.5c`$. The shock compression ratio is $`R=5.11`$. The particle spectra are measured at the shock for three different magnetic field perturbation amplitudes and the flat (Fig. 1a) or the Kolmogorov (Fig. 1b) wave power spectrum. One can note that the particle spectral indices deviate from the small amplitude results of the pitch angle scattering model \[3-6\]. In addition, the increasing magnetic field perturbations can produce non-monotonic changes of the particle spectral index – the feature which has not been discussed for parallel shocks so far. Analogously to oblique shock waves \[cf. an accompanying paper in this volume\], our particle spectra are non power-law in the full energy range and the shape of the spectrum vary with the amplitude of turbulence and the wave power index.

The non-monotonic variation of the spectral index with the turbulence amplitude results from modifications of the particle acceleration process at the shock. The long-wave finite-amplitude perturbations produce locally oblique mag-

netic field configurations and lead to occurrence of particles reflected from the compressed field downstream of the shock. Probability of reflection depends on the turbulence amplitude and the amount of field perturbations with wavelengths larger than the resonance wavelength for a given particle, as presented in Figs. 1c and 1d. For $`\delta B/B_{0,1}=1.0`$ the reflection probability is higher as compared to the other perturbation amplitudes considered and the particle spectrum is flatter. For the chosen by us smaller ($`\delta B/B_{0,1}=0.3`$) and larger ($`\delta B/B_{0,1}=3.0`$) turbulence amplitudes the reflection and transmission probability do not differ considerably, which results in the similar values of the spectral indices. One can note in this place, that the obtained spectra for the Kolmogorov case seem to exhibit a continuous slow change of inclinations. Thus the fitted power-laws depend to some extent on the energy range chosen for the fit. One can also note a steep part of the spectrum at low energies for $`\delta B/B_{0,1}=0.3`$ in Fig. 1b.

The presented reflection (transmission) probabilities decrease (increase) at high particle energies due to a limited dynamic range of the magnetic field turbulence. The locally oblique field configurations are mainly formed by long-wave perturbations ($`k<k_{res}`$) \[cf. 9\]. For high energy particles with $`k_{res}<k_{min}`$ there are no respectively long waves and the upstream particles can be only transmitted downstream of the shock. In these conditions the acceleration process would converge to the ‘classic’ parallel shock acceleration model, but in our simulations particles move far to the introduced escape boundary forming a cut-off.

### 4. Summary

The simulations of the first-order Fermi acceleration process acting at parallel relativistic shock waves are presented. In the presence of finite-amplitude perturbations the particle spectral indices vary in a non-monotonic way with the turbulence amplitude due to non trivial character of the particle interaction with the shock, including particle reflections. This feature is also expected to lead to some decrease of the particle acceleration timescale in comparison to previous estimates \[7, cf. 8\].

The work was supported by the Polish State Committee for Scientific Research in 2002-2004 as a research project 2 P03D 008 23 (JN) and in 2002-2005 as a solicitated research project PBZ-KBN-054/P03/2001 (M0).

### References

1. Bednarz J., Ostrowski M. 1996, MNRAS 283, 447
2. Bednarz J., Ostrowski M. 1998, Phys. Rev.Lett. 80, 3911
3. Heavens A., Drury L’O.C. 1988, MNRAS 235, 997
4. Kirk J.G., Heavens A. 1989, MNRAS 239, 995
5. Kirk J.G., Schneider P. 1987, ApJ 315, 425
6. Kirk J.G., Schneider P. 1987, ApJ 322, 256
7. Lagage P.0., Cesarsky C.J. 1983, A&Ap 118, 223
8. Ostrowski M. 1988, MNRAS 233, 257
9. Ostrowski M. 1988, A&Ap 206, 169
10. Ostrowski M. 1993, MNRAS 264, 248