text:
# Multiple scattering of polarized radiation by non-spherical grains: first results

## 1 Introduction

It is now well established that the polarization of optical and near-infrared radiation from young and evolved stellar objects, reflection nebulae, and active galactic nuclei is mainly caused by dust grains. In many cases the observed polarization can be satisfactorily interpreted by light scattering on spherical grains. In particular, the orientation of the polarization vectors is often used for the investigation of the dust distribution and for identifying the location of the embedded illuminating source(s).

However, even in the case of simple objects like ordinary reflection nebulae, deviations of polarization vectors from the direction perpendicular to a star were discovered more than three decades ago (Elvius & Hall eh67 (1967)). These deviations are displayed in the outer filamentary parts of the Merope nebula and may be easily explained by light scattering on non-spherical grains aligned by a magnetic field along the filaments. Non-centrosymmetric polarization patterns have been observed in bipolar and cometary nebulae (Scarrott et al. sdw89 (1989)), young stellar objects (Hajjar & Bastien hb96 (1996)), evolved stars (Kastner & Weintraub kw96 (1996)), and are also clearly seen in polarization maps of comets (Dollfus & Suchail ds87 (1987)). They may be attributed to non-spherical grains, although they could be caused by multiple scattering on spherical particles as well. Another effect which may be related to the light scattering by non-spherical grains is the wavelength dependence of the positional angle of polarization observed in red giants, AGB stars, and bipolar reflection nebulae (see, e.g., Johnson & Jones jj91 (1991)). The variations from blue to red may reach $`20\mathrm{°}60\mathrm{°}`$ and it is very difficult or even impossible to interpret this behaviour using spherical grains only.

Recently, very high degrees of circular polarization of scattered light in the Orion molecular cloud were measured by Chrysostomou et al. (cgm00 (2000)). The authors suggest that the circular polarization is produced by aligned non-spherical grains. Multiple scattering of radiation by spherical or randomly oriented non-spherical grains results in a much smaller circular polarization degree than observed. In addition, the interstellar polarization and polarized thermal emission phenomena prove that non-spherical grains exist in the interstellar medium. These effects arise because of dichroic extinction/emission of radiation by aligned non-spherical grains and were modeled with spheroidal grains (Voshchinnikov v90 (1990); Kim & Martin km95 (1995); Onaka o00 (2000)).

Up to now, spheroidal grains have been used for the interpretation of scattered radiation only in the case of single scattering (Voshchinnikov v98 (1998); Gledhill & McCall gm00 (2000)). The numerical simulation of the polarized radiation transfer through a medium with non-spherical grains, including multiple scattering, is extremely difficult. In the case of spherical grains nearly all previous simulations are based on Monte Carlo simulations. This method was used for the interpretation of polarimetric observations of various stars (see, e.g., Daniel d80 (1980); Voshchinnikov & Karjukin vk94 (1994)) as well as the production of polarimetric maps of different extended objects (e.g., Bastien & Ménard bm88 (1988); Fischer et al. fhy94 (1994); Wolf et al. wfp98 (1998)). The most recent, major achievement of the Monte Carlo technique — in respect to the solution to the radiative transfer (RT) problem — is the development of RT codes that allow to calculate the spectral energy distribution of multi-dimensional dust configurations self-consistently (Wolf et al. whs99 (1999); Wolf & Henning wh00 (2000)).

The main aim of this paper is to provide the first radiative transfer simulations including scattering, extinction, absorption, and re-emission of radiation by aligned spheroidal dust grains. Our formalism will be applied to several simple model configurations. Because of the complexity of the problem, we decided to consider the separate mechanisms step by step. In this paper, the basic theory is presented and the main numerical features of the simulation are described (Sect. 2). Furthermore, we consider multiple scattering of light by spheroidal grains and the resulting linear and circular polarization (Sect. 3 and 4).

## 2 Basics

### 2.1 Dust grain parameters

Let us suppose that the dust grains are prolate or oblate homogeneous spheroids with an aspect ratio $`a/b`$, where the quantities $`a`$ and $`b`$ are the major and minor semiaxes of a spheroid, respectively. The particle size is characterized by the radius $`r_\mathrm{V}`$ of a sphere with the same volume as that of the spheroid (equivolume radius). The major semiaxis $`a`$ of the spheroid is connected with $`r_\mathrm{V}`$ by

$$a=r_\mathrm{V}\left(\frac{a}{b}\right)^{2/3}$$
(1)
for prolate spheroids and

$$a=r_\mathrm{V}\left(\frac{a}{b}\right)^{1/3}$$
(2)
for oblate spheroids. The optical properties of a particle also depend on the complex refractive index of the grain material $`m_\lambda `$ and the angle $`\alpha `$ between the rotation axis of the spheroid and the wave-vector direction (direction of the incident radiation; $`0\mathrm{°}\alpha 90\mathrm{°}`$, see Fig. 1).

### 2.2 Radiative transfer concept

In this section we outline the main ideas of our solution of the RT problem taking into account multiple scattering and dichroic absorption, and reemission by spheroidal grains. Our solution is based on the Monte Carlo method which is described only briefly here. For a more general overview about RT based on the Monte Carlo method, we refer to the papers of Wolf et al. (whs99 (1999)) and Wolf & Henning (wh00 (2000)) (see also Fischer et al. fhy94 (1994)).

In our solution to the RT problem, the radiation energy is partitioned into so-called weighted photons (in the following we will use the term “photon”). The energy, intensity, and polarization of a weighted photon are described by its wavelength $`\lambda `$ and the Stokes vector $`\widehat{I}=(I,Q,U,V)^\mathrm{T}`$. The transformation of the Stokes vector by the $`i`$th scattering can be described by a $`4\times 4`$ (Müller) matrix $`\widehat{F}`$:

$$\left(\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right)_i\widehat{F}\left(\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right)_{i1}.$$
(3)
The main features distinguishing the light scattering by aligned non-spherical particles and spherical particles are:

1. Azimuthal dependence of the scattered radiation;
2. Linear polarization in the forward and backward directions;
3. Deviation of the positional angle of linear polarization after first scattering from the direction perpendicular to the illuminating source and
4. Circular polarization after first scattering.

These features result from the fact that in the general case of non-spherical particles — in contrast to spherical grains — all 16 elements of the scattering matrix are non-zero (see Voshchinnikov & Farafonov vf93 (1993)).

The spatial distribution of the illuminating star(s) and the dust grains is defined inside a convex model space. The free path length $`l`$ of a photon (distance between the point of emission and the first point of scattering/absorption or two points of scattering/absorption) is determined by the dust density distribution $`n_\mathrm{d}(R_i)`$ at the coordinate $`R_i`$ and the extinction cross-section for the transmitted radiation $`\stackrel{~}{C}_{\mathrm{ext}}(R_i)`$

$$\stackrel{~}{\tau }_{\mathrm{ext}}=\underset{i=1}{\overset{i_{\mathrm{end}}}{}}n_\mathrm{d}(R_i)\stackrel{~}{C}_{\mathrm{ext}}(R_i)\mathrm{\Delta }l_i.$$
(4)
The quantity $`i_{\mathrm{end}}`$ represents the number of integration steps along the free path length, and the quantity $`\mathrm{\Delta }l_i`$ is the corresponding geometrical step width. The free mean path length

$$l=\underset{i=1}{\overset{i_{\mathrm{end}}}{}}\mathrm{\Delta }l_i$$
(5)
can be derived from the optical depth $`\stackrel{~}{\tau }_{\mathrm{ext}}`$ as follows

$$\stackrel{~}{\tau }_{\mathrm{ext}}=\mathrm{ln}(1\zeta ).$$
(6)
Here, $`\zeta `$ is a random number uniformly distributed on the interval $`[0,1]`$ (see Wolf et al. whs99 (1999) for details).

The extinction cross-section $`\stackrel{~}{C}_{\mathrm{ext}}`$ is calculated taking into account the state of polarization of the incident radiation (see, e.g., Martin m74 (1974)):

$`\stackrel{~}{C}_{\mathrm{ext}}={\displaystyle \frac{1}{2}}(C_{\mathrm{ext}}^{\mathrm{TM}}+C_{\mathrm{ext}}^{\mathrm{TE}})+{\displaystyle \frac{1}{2}}(C_{\mathrm{ext}}^{\mathrm{TM}}C_{\mathrm{ext}}^{\mathrm{TE}})\mathrm{cos}2\mathrm{\Psi }{\displaystyle \frac{Q_0}{I_0}}+`$
$`{\displaystyle \frac{1}{2}}(C_{\mathrm{ext}}^{\mathrm{TM}}C_{\mathrm{ext}}^{\mathrm{TE}})\mathrm{sin}2\mathrm{\Psi }{\displaystyle \frac{U_0}{I_0}}.`$ (7)
Here, $`\mathrm{\Psi }`$ is the angle between the particle frame and the laboratory frame. $`(I_0,Q_0,U_0,V_0)`$ is the Stokes vector of the incident light. The extinction cross-sections $`C_{\mathrm{ext}}^{\mathrm{TE},\mathrm{TM}}`$ are connected with the efficiency factors $`Q_{\mathrm{ext}}^{\mathrm{TE},\mathrm{TM}}`$ by

$$C_{\mathrm{ext}}^{\mathrm{TE},\mathrm{TM}}=GQ_{\mathrm{ext}}^{\mathrm{TE},\mathrm{TM}},$$
(8)
where the quantity $`G`$ is the geometrical cross-section of the dust grain (the area of the particle shadow). The superscripts TM and TE are related to the two cases of the polarization of the incident radiation (TM and TE modes). The geometrical cross-section $`G`$ can be derived as follows

$$G(\alpha )=\pi r_\mathrm{V}^2\left(\frac{a}{b}\right)^{2/3}\left[\left(\frac{a}{b}\right)^2\mathrm{sin}^2\alpha +\mathrm{cos}^2\alpha \right]^{1/2}$$
(9)
for a prolate spheroid and

$$G(\alpha )=\pi r_\mathrm{V}^2\left(\frac{a}{b}\right)^{2/3}\left[\left(\frac{a}{b}\right)^2\mathrm{sin}^2\alpha +\mathrm{cos}^2\alpha \right]^{1/2}.$$
(10)
for an oblate spheroid. The efficiency factors $`Q_{\mathrm{ext}}^{\mathrm{TM},\mathrm{TE}}`$ are calculated from the solution to the light scattering problem for spheroids obtained by the method of separation of variables (see Voshchinnikov & Farafonov vf93 (1993) for details).

We would like to note that the extinction cross-section $`\stackrel{~}{C}_{\mathrm{ext}}`$ (see Eq. (2.2)) has to be determined at each integration step along the photon path (see Eq. (4)). In the case of the RT in a medium with aligned non-spherical grains, the step width along the photon free path is determined not only by the scale of the density structure (as for spherical particles) but also by the degree of spatial variations of the grain alignment. Owing to this and the more complex treatment of the scattering process compared with spherical grains, the total computation run-time is by an order of magnitude larger for aligned spheroids than for spheres.

Finally, depending on the model geometry and the dust density distribution, the photon leaves the model space after a certain number of interactions with dust grains (scattering, absorption). Then, it will be “observed” by array-detector-like planes arranged around the model space. The Stokes vector has to be projected onto the respective plane. For a detailed description of the projection formalism we refer to Fischer et al. (fhy94 (1994)).

### 2.3 Scattering

The solution to the problem of light scattering by spheroidal particles is given in a particle reference system related to the rotation axis of a spheroid (in spherical coordinates $`\theta ,\phi `$), while the RT is considered in the laboratory reference system (in spherical coordinates $`\mathrm{\Theta },\mathrm{\Phi }`$; see Fig. 1). The probability of a photon of being scattered into the direction ($`\mathrm{\Theta },\mathrm{\Phi }`$) is proportional to the value of the first element of the scattering matrix $`F_{11}(\alpha ,\mathrm{\Theta },\mathrm{\Phi })`$. We transformed the scattering matrix from the particle to the laboratory frame. Then, the scattering matrix is calculated using the standard expressions (Bohren & Huffman bh83 (1983)). In order to check the transition from the particle system to the laboratory one, the transformation was also done for light scattering by spheroidal particles in the quasistatic approximation where the solution can be written in both coordinate systems (see Farafonov et al. fip01 (2001) for details).

The elements of the scattering matrix were normalized using the standard procedure (Bohren & Huffman bh83 (1983))

$`{\displaystyle \frac{2}{k^2(C_{\mathrm{sca}}^{\mathrm{TM}}+C_{\mathrm{sca}}^{\mathrm{TE}})}}{\displaystyle _0^{2\pi }}{\displaystyle _0^\pi }F_{11}(\alpha ,\mathrm{\Theta },\mathrm{\Phi })\mathrm{sin}\mathrm{\Theta }\mathrm{d}\mathrm{\Theta }\mathrm{d}\mathrm{\Phi }=1,`$ (11)
where $`C_{\mathrm{sca}}`$ is the scattering cross-section and $`k=2\pi /\lambda `$ the wavenumber. The correctness of these computations was verified using the relationships for matrices, describing scattering by small particles (see Hovenier & van der Mee hv00 (2000) and references therein).

The calculation of the elements of the scattering matrix $`\widehat{F}(m,r_\mathrm{V},a/b,\alpha ,\mathrm{\Theta },\mathrm{\Phi })`$ and efficiency factors $`Q_{\mathrm{sca},\mathrm{ext}}^{\mathrm{TM},\mathrm{TE}}(m,r_\mathrm{V},a/b,\alpha )`$ for a single particle with the fixed equivolume radius $`r_V`$, the aspect ratio $`a/b`$, and the refractive index $`m`$ requires several minutes of computer time111A Pentium II 450 MHz processor was used for our calculations.. These quantities were precalculated before starting the RT modelling. In order to guarantee that the main features (local minima/maxima) of the scattering matrix elements (and/or derived quantities, for instance, the polarization of singly scattered light) are present in the discrete distributions as a function of $`\alpha `$, $`\mathrm{\Theta }`$, and $`\mathrm{\Phi }`$ (see Figs. 2 and 3), we use step widths of $`\mathrm{\Delta }\alpha =\mathrm{\Delta }\mathrm{\Theta }=\mathrm{\Delta }\mathrm{\Phi }=2\mathrm{°}`$ in our simulations (see Fig. 3).

### 2.4 Absorption and Reemission

The absorption is described by the albedo $`\mathrm{\Lambda }`$ of the particle which is defined as the ratio of the scattering and extinction cross sections. It is equal to

$$\mathrm{\Lambda }^{\mathrm{TM},\mathrm{TE}}(m,r_\mathrm{V},a/b,\alpha )=\frac{C_{\mathrm{sca}}^{\mathrm{TM},\mathrm{TE}}(m,r_\mathrm{V},a/b,\alpha )}{C_{\mathrm{ext}}^{\mathrm{TM},\mathrm{TE}}(m,r_\mathrm{V},a/b,\alpha )}$$
(12)
for totally polarized light and

$$\mathrm{\Lambda }^{(0)}(m,r_\mathrm{V},a/b,\alpha )=\frac{C_{\mathrm{sca}}^{\mathrm{TM}}+C_{\mathrm{sca}}^{\mathrm{TE}}}{C_{\mathrm{ext}}^{\mathrm{TM}}+C_{\mathrm{ext}}^{\mathrm{TE}}}$$
(13)
for non-polarized light.

While in the case of absorption by spherical particles, the Stokes vector $`\widehat{I}`$ can be multiplied by a scalar value of the albedo, for spheroidal particles the following $`4\times 4`$ matrix is needed (in the particle frame):

$$\widehat{\mathrm{\Lambda }}_{\mathrm{part}.}=\left(\begin{array}{cccc}\stackrel{~}{l}_1& \stackrel{~}{l}_2& 0& 0\\ \stackrel{~}{l}_2& \stackrel{~}{l}_1& 0& 0\\ 0& 0& \stackrel{~}{l}_1& 0\\ 0& 0& 0& \stackrel{~}{l}_1\end{array}\right)_{\mathrm{part}.}.$$
(14)
Here, the quantities $`\stackrel{~}{l}_1`$ and $`\stackrel{~}{l}_2`$ are defined as

$`\stackrel{~}{l}_1=\mathrm{\Lambda }={\displaystyle \frac{\mathrm{\Lambda }^{\mathrm{TM}}}{1+C_{\mathrm{ext}}^{\mathrm{TM}}/C_{\mathrm{ext}}^{\mathrm{TE}}}}+{\displaystyle \frac{\mathrm{\Lambda }^{\mathrm{TE}}}{1+C_{\mathrm{ext}}^{\mathrm{TE}}/C_{\mathrm{ext}}^{\mathrm{TM}}}},`$
$`\stackrel{~}{l}_2={\displaystyle \frac{\mathrm{\Lambda }^{\mathrm{TM}}}{1+C_{\mathrm{ext}}^{\mathrm{TM}}/C_{\mathrm{ext}}^{\mathrm{TE}}}}{\displaystyle \frac{\mathrm{\Lambda }^{\mathrm{TE}}}{1+C_{\mathrm{ext}}^{\mathrm{TE}}/C_{\mathrm{ext}}^{\mathrm{TM}}}}.`$ (15)
Before applying the albedo matrix $`\widehat{\mathrm{\Lambda }}_{\mathrm{part}.}`$, it has to be rotated into the particle frame using the rotation matrix $`\widehat{R}`$

$$\widehat{R}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}2\mathrm{\Psi }& \mathrm{sin}2\mathrm{\Psi }& 0\\ 0& \mathrm{sin}2\mathrm{\Psi }& \mathrm{cos}2\mathrm{\Psi }& 0\\ 0& 0& 0& 1\end{array}\right).$$
(16)
The absorption process can then be described by the expression

$$\widehat{I}_{\mathrm{after}}=\widehat{R}^1\widehat{\mathrm{\Lambda }}_{\mathrm{part}.}\widehat{R}\widehat{I}_{\mathrm{before}},$$
(17)
whereby $`\widehat{I}_{\mathrm{before}}`$ ($`\widehat{I}_{\mathrm{after}}`$) is the Stokes vector of the radiation before (after) absorption.

In the case of (re)emission of radiation by a spheroidal particle, the Stokes vector of the photon can be written as follows (particle frame):

$$\left(\begin{array}{c}I_0\\ Q_0\\ U_0\\ V_0\end{array}\right)_{\mathrm{part}.}\left(\begin{array}{c}(C_{\mathrm{abs}}^{\mathrm{TM}}+C_{\mathrm{abs}}^{\mathrm{TE}})/2\\ (C_{\mathrm{abs}}^{\mathrm{TM}}C_{\mathrm{abs}}^{\mathrm{TE}})/2\\ 0\\ 0\end{array}\right)_{\mathrm{part}.}B_\lambda (T_\mathrm{d}).$$
(18)
Here, $`C_{\mathrm{abs}}`$ is the absorption cross-section, $`B_\lambda (T_\mathrm{d})`$ the Planck function and $`T_\mathrm{d}`$ the grain temperature.

## 3 Test simulations

The best way to test any RT code is the comparison of its results with those obtained with another RT code/technique or an analytical solution. Because the RT code presented here is the first that provides the treatment of multiple scattering of polarized light by aligned non-spherical particles, a lack of suitable test cases for comparison is evident. Nevertheless, in the following we present test calculations which confirm the correct implementation of the new scattering mechanism into our RT code.

For test simulations we used particles with the following parameters:

* shape of the grains: oblate,
* refractive index: $`m=1.7+0.03i`$,
* equivolume radius: $`r_\mathrm{V}=0.1\mu `$m,
* aspect ratio: $`a/b=2.0`$,
* wavelength of incident radiation: $`\lambda =0.628\mu `$m.

Then the grain size parameter is $`x_\mathrm{V}=2\pi r_\mathrm{V}/\lambda =1`$.

### 3.1 Single scattering in a geometrically and optically thin hemisphere

Firstly, we consider single scattering in a thin layer of oblate spheroids at the surface of an hemisphere (see Fig. 4).

The ratio of the inner radius ($`r_{\mathrm{in}}`$) of the shell to the outer one ($`r_{\mathrm{out}}`$) amounts to $`r_{\mathrm{in}}/r_{\mathrm{out}}=0.99`$. The non-rotating particles show a static picket fence orientation with their rotation axes being aligned parallel to the $`z`$-axis of the global coordinate system of the model space. The density profile $`n(r)`$ in the shell follows a power law:

$$n(r)r^1.$$
(19)
The optical thickness of the shell at the wavelength of the simulated RT ($`\lambda =0.628\mu `$m) is 0.2 along the $`z`$-axis, i.e. $`\tau _\mathrm{z}(0.628\mu \mathrm{m})=0.2`$. Because of the chosen grain alignment, $`\tau _\mathrm{z}`$ is proportional to the extinction cross-section in the case of light incident parallel to the grain rotation axis, $`\tau _0=\tau _\mathrm{z}C_{\mathrm{ext}}(0\mathrm{°})`$. In other directions, the optical thickness of the shell may be found as $`\tau _\alpha =\tau _\mathrm{z}C_{\mathrm{ext}}(\alpha )/C_{\mathrm{ext}}(0\mathrm{°})`$ and, in particular, along the line of sight ($`x`$-axis) and the $`y`$-axis we have $`\tau _\mathrm{x}=\tau _\mathrm{y}=\tau _{90\mathrm{°}}=\tau _\mathrm{z}C_{\mathrm{ext}}(90\mathrm{°})/C_{\mathrm{ext}}(0\mathrm{°})`$. Using Eq. 2.2 for non-polarized stellar radiation and the standard model described above, we have $`\tau _{90\mathrm{°}}0.35`$ (prolate grains) and $`\tau _{90\mathrm{°}}0.13`$ (oblate grains), if $`\tau _0(0.628\mu \mathrm{m})=0.2`$. For more elongated or flattened particles ($`a/b=10`$), these values become equal to $`\tau _{90\mathrm{°}}0.68`$ and 0.082 for prolate and oblate grains, respectively.

In contrast to polarization maps obtained for model configurations with the same geometry, but containing spherical dust grains, the polarization pattern is not centro-symmetric (see Fig. 5\[A\]). This is because of the azimuthal asymmetry of light scattering by aligned non-spherical particles. The net polarization in this case is not equal to zero. We compare our results with those obtained by Matsumura & Seki (ms96 (1996))222The numerical data in tabular form were provided by Matsumura (m99 (1999))..

In Fig. 5 the comparison of the polarization degree as a function of the projected distance from the centre along three different lines in the polarization map is shown. These graphs as well as the polarization pattern (see for comparison Fig. 2 from Matsumura & Seki ms96 (1996)) are in perfect agreement.

Qualitatively, the polarization pattern is quite similar to that of a young stellar object surrounded by a spherical shell and an optically thick disk seen edge-on. There the parallel polarization pattern is due to light scattering on the disk surface while the centro-symmetric polarization pattern is due to light scattering in the circumstellar environment. Corresponding RT simulations with spherical dust grains have been performed by M$`\stackrel{´}{\mathrm{e}}`$nard (m88 (1988)) and Fischer et al. (fhy94 (1994)). The most interesting feature of the pattern is the presence of polarization null points (see Figs. 5\[A\] and 5\[B1\]). They have been observed in a number of objects (see discussion in Sect. 4.2.3). Sometimes, these null points are interpreted as the indicators for the “outer edge” of a circumstellar disk (transition from the optically thick disk to the optically thin shell; see Fischer et al. fhy94 (1994)). Although this explanation is very plausible, the simulations presented here show that uniformly aligned grains (possibly by an external magnetic field) may produce a similar polarization pattern. A detailed discussion of the polarization null points as a function of the model parameters is given in Sect. 4.

### 3.2 Multiple scattering

To test the RT with multiple scattering, we consider model configurations for which the dust density distribution, the dust grain orientation, and the irradiation characteristic are centro-symmetric with respect to the observer’s position. The resulting polarization pattern is therefore also centro-symmetric: the polarization vectors are oriented perpendicular or parallel to the radius vector of the configuration, in dependence of the refractive index and the grain size parameter $`x_\mathrm{V}`$. Thus, the integral (net) linear and circular polarization have to be zero. Because this property of centro-symmetric configurations does not depend on its optical thickness, it can be used as a test of the single and multiple scattering case.

We considered the following dust configurations containing either oblate or prolate grains:

1. Spherical dust shell with a point-like radiation source in the centre and particle rotation axes either parallel or perpendicular to the radius vector;
2. Cylindrical dust shell with a point-like radiation source on its rotation axis and particle rotation axes either parallel or perpendicular to the radius vector which is oriented perpendicular to the rotation axis of the cylinder.

While for the first configuration the net polarization is equal to zero for every position of the observer, for the second configuration this is the case only if the observer’s line of sight is oriented perpendicular to one of the two cylinder footprints. In all cases we arrived at the expected results.

## 4 Results and discussion

In Sect. 4.1, light scattering in a geometrically thin hemisphere and sphere and a geometrically thick spherical shell is described. This preliminary modelling provides a “smooth” transition from the test models discussed in Sect. 3 to more realistic dust configurations. The basic modelling is presented in Sect. 4.2, where a geometrically thick shell is considered and the resulting polarization as a function of the main model parameters is investigated.

If one considers the RT in a medium consisting of spherical dust grains, scattering is the only mechanism to produce or to modify the polarization of the incident radiation. In the case of aligned spheroidal grains, the processes of photon absorption and re-emission may also fabricate or change polarization. Because of the obvious complexity of the problem, we separated the consideration of these three mechanisms of the polarization formation. The aim of the simulations discussed below was to outline the basic features of polarization arising from the light scattering by spheroidal grains. The additional influence of absorption and re-emission will be discussed in future publications.

For this reason, we apply the following approximations:

1. The grain albedo is assumed to be that of spherical grains with the equivolume radius $`r_\mathrm{V}`$ and is therefore a scalar. Then, the absorption process decreases the intensity of the incident radiation but does not modify the state of polarization.
2. We do not consider thermal re-emission of radiation.

If not stated differently, for the simulation presented in this section the dust grain and model parameters as described in Sect. 3 have been used. The observer’s position is inclined by $`i=90^\mathrm{o}`$ to the $`z`$-axis ($`\tau _\mathrm{z}(0.628\mu \mathrm{m})=0.2`$) of the global coordinate system of the model space (see Fig. 4).

### 4.1 Simple models

Let us consider again the geometrically thin hemispherical model. While the polarization pattern in Fig. 5 is shown only for a position perpendicular to the $`z`$-axis ($`i=90\mathrm{°}`$), Fig. 6 demonstrates how the behaviour of the degree of linear polarization in the midplane of the polarization map appears for different inclinations of the $`z`$-axis to the line of sight. It turns out that the local polarization minimum (the “polarization null point” or the reversal of polarization in the case $`i=90^\mathrm{o}`$) are present at large inclination angles $`i`$ only. It is important to note that the polarization pattern (and therefore also the net polarization degree) for aligned non-spherical particles depends on the viewing angle even in the case of a one-dimensional density distribution. In contrast to this, no dependence of the polarization on the inclination angle exists for spherical particles.

As a first step away from the very artificial configuration described above, we consider the RT in a thin spherical shell instead of the hemisphere shown in Fig. 4. Then, backward scattering in the hemisphere opposite to the observer’s position (where $`\mathrm{\Theta }=90\mathrm{°}+\alpha `$) has to be taken into account. The results are shown in Fig. 7. The addition of the backscattered radiation shifts the position of the minimum of the linear polarization to a smaller projected radius.

In Fig. 8, the circular polarization for this model of a geometrically and optically thin dust shell is shown. The circular polarization is equal to zero in the midplane and along the line through the centre perpendicular to the midplane. The net circular polarization from the configuration is also equal to zero. A counterpart of the linear polarization null points in the midplane could not be found. The maximum circular polarization is observed at the ends of lines inclined to the $`z`$-axis by $`\pm 45\mathrm{°}`$. These points correspond to the light scattering by particles with values of $`\alpha =45\mathrm{°}`$, $`\mathrm{\Theta }=90\mathrm{°}`$ and $`\mathrm{\Phi }=90\mathrm{°}(270\mathrm{°})`$.

It turned out that the pattern of the linear polarization strongly depends on the inner radius of the spherical shell. This is demonstrated in Fig. 9, where the degree of linear polarization in the midplane is shown for three ratios of the inner to the outer radius of the shell. Decreasing this ratio (and therefore decreasing the relative size of the inner cavity around the illuminating star), the characteristic polarization null points shift towards the star. It has to be noted, that generally the polarization null point is not located at the position of the projected inner radius. The location, relative depth and width of the area where the polarization changes the sign depends on the geometrical thickness of the shell and grain parameters (see Sect. 4.2.2). The power-law density profile used (see Eq. (19)) leads to a preferred photon scattering near a narrow region at the inner boundary of the shell. The shift of the polarization null point is therefore mainly caused by the change of the inner radius $`r_{\mathrm{in}}`$.

### 4.2 Envelope models

Now we will consider geometrically thick spherical shells. To represent such models, the ratio of the inner radius of the shell to the outer radius is fixed at $`r_{\mathrm{in}}/r_{\mathrm{out}}=0.5`$.

#### 4.2.1 Transition from optically thin to optically thick shell

In the case of a small optical thickness, the main contribution to the resulting polarization is from first or second scattering events (Voshchinnikov & Karjukin vk94 (1994); Voshchinnikov et al. vgk95 (1995)). As one can see from Fig. 10, the polarization $`P_\mathrm{l}(r_y/r_{\mathrm{out}})`$ in the midplane decreases with increasing optical thickness. The same behaviour was found for the polarization perpendicular to the midplane $`P_\mathrm{l}(r_z/r_{\mathrm{out}})`$. In general, the position of the local polarization minimum is not strongly influenced by the optical thickness.

Integrating the images for the Stokes parameters $`I,Q,U`$, we estimate the total linear polarization from our spherical dust configuration as one observes from the position perpendicular to the $`z`$-axis. Because of the symmetry, it is the same along the $`x`$ and $`y`$-axes

$$P_{\mathrm{l},\mathrm{\hspace{0.17em}90}\mathrm{°}}=\frac{F_{\mathrm{sca},\mathrm{\hspace{0.17em}90}\mathrm{°}}^\mathrm{p}}{F_{}\mathrm{exp}\left[\tau _{90\mathrm{°}}^{\mathrm{ext}}\right]+F_{\mathrm{sca},\mathrm{\hspace{0.17em}90}\mathrm{°}}},$$
(20)
where $`F_{\mathrm{sca},\mathrm{\hspace{0.17em}90}\mathrm{°}}`$ and $`F_{\mathrm{sca},\mathrm{\hspace{0.17em}90}\mathrm{°}}^\mathrm{p}`$ are the scattered and polarized scattered radiation from the configuration. For the models presented in Fig. 10, the resulting polarization is $`P_\mathrm{l}=4\%`$, 9%, and 8%, if $`\tau _\mathrm{z}=`$ 0.2, 1 and 2, correspondingly. Excluding the stellar radiation (first term in the denominator in Eq. (20)), we can find the maximum degree of the observed linear polarization

$$P_{\mathrm{l},\mathrm{\hspace{0.17em}90}\mathrm{°},\mathrm{max}}=\frac{F_{\mathrm{sca},\mathrm{\hspace{0.17em}90}\mathrm{°}}^\mathrm{p}}{F_{\mathrm{sca},\mathrm{\hspace{0.17em}90}\mathrm{°}}}.$$
(21)
It is equal to $`P_{\mathrm{l},\mathrm{\hspace{0.17em}90}\mathrm{°},\mathrm{max}}=25\pm 1\%`$, $`15\pm 1\%`$, and $`9.1\pm 1\%`$, if $`\tau _\mathrm{z}=`$ 0.2, 1 and 2, correspondingly.

#### 4.2.2 Dependence on dust grain parameters

Figure 11 shows the midplane polarization $`P_\mathrm{l}(r_\mathrm{y}/r_{\mathrm{out}})`$ for oblate and prolate grains of three different equivolume radii ($`r_\mathrm{V}`$ = 0.05$`\mu `$m, 0.1$`\mu `$m, 0.2$`\mu `$m). The behaviour of polarization for particles with $`r_\mathrm{V}=0.05\mu \mathrm{m}`$ and $`0.1\mu \mathrm{m}`$ is rather similar, while the polarization produced by the largest particles ($`r_\mathrm{V}=0.2\mu `$m) is the smallest and is almost independent of the projected radius. This feature reflects the fact of decreasing polarizing efficiency of large particles (Voshchinnikov et al. vihf00 (2000)). Because the manner of polarization for prolate grains in the direction perpendicular to the midplane is similar to that for oblate grains in the midplane and vice versa (see also Fig. 13), the resulting polarization for particles of both types is quite similar. The net polarization — excluding the direct, non-scattered stellar radiation — amounts to $`P_{\mathrm{l},\mathrm{\hspace{0.17em}90}\mathrm{°},\mathrm{max}}=`$ 20 % (16 %), 25 % (20 %) and 11 % (10 %) for oblate (prolate) grains with equivolume radii $`r_\mathrm{V}=0.05\mu \mathrm{m}`$, $`0.1\mu \mathrm{m}`$ and $`0.2\mu \mathrm{m}`$, respectively.

The distribution of the linear polarization $`P_\mathrm{l}(r_\mathrm{y}/r_{\mathrm{out}})`$ in the midplane of the spherical shell is plotted in Fig. 12 for particles with four different aspect ratios $`a/b`$ ($`a/b`$ = 1.1, 2, 4, 10). The polarization minima in the midplane ($`r_\mathrm{z}=0`$) appear for oblate grains only. Their position shifts to outer projected radii $`r_\mathrm{y}/r_{\mathrm{out}}`$ with increasing aspect ratio. If $`a/b1`$, the location of the polarization null point leads to $`r_\mathrm{y}/r_{\mathrm{out}}0`$. This is consistent with the behaviour of the degree of linear polarization for spherical grains where the null polarization can be found at the centre for the scattering angle $`\mathrm{\Theta }=0\mathrm{°}`$. The resulting polarization for oblate and prolate grains of different shape is given in Table 1. Note that the polarization produced by oblate particles is larger in all cases.

The intensity maps of our configuration ($`i=90^\mathrm{o}`$) show an additional effect. In Fig. 13, the images with overlaid polarization pattern, corresponding to the intensity and linear polarization of scattered light, are shown for oblate and prolate grains with the aspect ratios $`a/b`$=1.1 and 2. In the case of nearly spherical grains ($`a/b=1.1`$), the ring-like structure is clearly seen for oblate as well as for prolate grains. Otherwise, in the case of $`a/b=2`$ the intensity maxima appear at the points of intersection with the $`z`$($`y`$)-axis for oblate (prolate) grains. In both cases, the corresponding axis is parallel to the minor ($`b`$) axis of grains. Increasing the aspect ratio $`a/b`$, this effect is strengthened. Note that the net polarization from the spherical configuration arises in the inner parts of the shell because the polarization pattern from the outer regions is centro-symmetric and the polarization would be cancelled. The direction of the resulting net polarization is parallel to the major ($`a`$) axis of grains that is similar to the case of dichroic polarization.

#### 4.2.3 Null points on polarization maps as indicators of scattering by non-spherical grains

The model discussed in this paper is still too simple in order to make a detailed comparison with observations. However, we will discuss qualitatively the behaviour of null points produced by scattering on non-spherical grains.

Polarization null points were observed in a variety of objects ranging from reflection nebula and circumstellar structures to external galaxies (e.g., Scarrott et al. sdw89 (1989); Scarrott et al. srs90 (1990); Gledhill & Scarrott gs89 (1989); Asselin et al. assetal96 (1996); Chrysostomou et al. chryetal96 (1996); Hajjar & Bastien hb96 (1996); Kastner & Weintraub kw96 (1996); Wood & Jones wj97 (1997); Gledhill et al. getal01 (2001)). They appear approximately symmetrically on either side of the central illuminating source333It should be mentioned that the measurable polarization — especially near the position of polarization null points — sensitively depends on the spatial resolution of the polarization maps. Therefore, it is more likely to find the polarization minima but not real null points.. In all cases the origin of null points was attributed to the combined action of two polarization mechanisms, where the first (initial or internal) mechanism produced linearly polarized radiation whereas the second (external) mechanism led to the compensation of this polarization. Three possible models were qualitatively discussed by Scarrott et al. (sdw89 (1989)):

(i) Polarization due to scattering by spherical grains and its cancellation via dichroic extinction in an external medium consisting of aligned non-spherical particles;

(ii) Polarization due to scattering by spherical grains in an optically thick internal medium and its cancellation after the next scattering by spherical grains in an external optically thin medium;

(iii) Polarization of light from the central source due to an unspecified mechanism and its compensation due to the second scattering by spherical grains in an external medium.

These models could be called “hybrid models” because they require two steps to produce the observed polarization null points. The model (i) is rather naturally applied to the modelling of polarization maps of galaxies (Wood & Jones wj97 (1997)). The multiple scattering model (ii) was used by Bastien & Ménard (bm88 (1988)) and Fischer et al. (fhy94 (1994)) for interpretation of polarization maps of young stellar objects. In this case, the centrosymmetric polarization pattern is obtained in outer optically thin parts of the shell, whereas the aligned polarization vectors appear in the dense inner parts of this shell. The polarized source model (iii) was discussed by Notni (n85 (1985)) and Gledhill (g91 (1991)). Using this model, it is possible to explain the polarization pattern observed in several bipolar and cometary nebulae (e.g., Gledhill g91 (1991)), but as it was shown by Clark et al. (cetal00 (2000)) that the models with spherical grains fail to produce even the 5 % circular polarization found in the reflection nebulosity around the Herbig Ae/Be star R CrA. In this case, the authors suggested scattering by aligned non-spherical grains as the operating mechanism.

In contrast to the hybrid models, different polarimetric effects can arise as a result of light scattering by aligned non-spherical particles. In our simple model, the location of polarization null points depends on the particle shape, size and the inner to outer radius ratio (see Figs. 9, 11 and 12). The wavelength dependence of the location of null points is not very strong and their position only slightly approaches the central (stellar) position when going from the V to K band as it follows from Fig. 14444In the case of polarization null points which have been observed at the location of the outer boundary of circumstellar disks around young stellar objects (transition from the optically thick disk to the optically thin shell), the same behaviour is expected since this boundary is slightly shifted inward with increasing wavelength..

A similar conclusion can be reached from the observed polarization maps of the variable Hubble’s nebula NGC 2261 (Scarrott et al. sdw89 (1989)) and the infrared reflection nebula GSS 30 (Chrysostomou et al. chryetal96 (1996)). The results of Asselin et al. (assetal96 (1996)) for V376 Cas show that the position of the null points is also shifted as a function of the binned box size used in averaging the polarization. However, it should be noted that near the null points the observed signal is usually very small and the errors are large.

Finally, we can indicate that the scattering by almost spherical grains is marked by a ring-like structure in the intensity maps whereas the null point are located inside this structure rather close to the stellar position (Fig. 13, left side). In the case of aligned non-spherical particles (Fig. 13, right side), the null points at the inner boundary of the spherical shell are accompanied by symmetric arc-like structures in the intensity maps. There is some similarity of such maps with those of small reflection nebulae published by Gledhill et al. (getal01 (2001)): the “spherical case” is seen in the source IRAS 19114+0002 while the “non-spherical case” occurs in the maps of IRAS 22223+4327, IRAS 19500-1709 and IRAS 17106-3046.

We should emphasize that in the presented model the origin of null points is caused by light scattering by non-spherical particles, but is not related to the cancellation of polarization as in hybrid models.

## 5 Conclusions

The paper contains the first solution to the RT problem including multiple light scattering in dust configurations with aligned non-spherical dust grains. Our main aim was to discuss the effects of light scattering by spheroidal dust grains in a simple (spherical) dust configuration, containing perfectly aligned particles.

The most remarkable features of the simulated linear polarization maps are the polarization null points where the reversal of polarization occurs. They appear when the grain alignment axis is perpendicular to the line of sight. Symmetrically to the polarization null points, we found maxima in the intensity maps.

In contrast to spherical grains, even single scattering by spheroidal grains may cause circular polarization. The maps of circular polarization have a sector-like structure with polarization maxima at the ends of a line inclined to the grain alignment axis by $`\pm 45\mathrm{°}`$.

Based on our theoretical investigations of light scattering, dichroic absorption, and (re)emission by spheroidal grains presented in Sect. 2, the next steps of the numerical radiative transfer simulations will include the consideration of light scattering by partly aligned rotating particles and the calculation of the polarized re-emission of grains in different dust configurations. This will immediately lead to an improved interpretation and therefore better understanding of mid-infrared and submillimeter polarization where re-emission and dichroic absorption by partly aligned non-spherical dust grains are the main polarization mechanisms. Thus, it provides a basis for investigations of the structure of the magnetic fields in these objects which are in most cases assumed to cause grain alignment.

###### Acknowledgements.

We appreciate the discussions of light scattering by non-spherical particles with Victor Farafonov and Bernhard Michel and wish to thank M. Matsumura for providing his results for comparison. We are grateful to the referee for constructive comments. This research was supported by the DFG grant Ste 605/10 within the program “Physics of star formation” and by grants of the Volkswagen Foundation and INTAS (Open Call 99/652).


format: text

text:
# Cosmic Ray Production of 6Li by Structure Formation Shocks in the Early Milky Way: A Fossil Record of Dissipative Processes during Galaxy Formation

## 1. Introduction

The light elements Li, Be and B are unique in that apart from 7Li, none can be synthesized appreciably in thermal environments such as stellar interiors, supernova envelopes or the standard Big Bang. Instead, the bulk of these elements are believed to arise from nonthermal nuclear reactions induced by cosmic rays (CRs). Their abundances observed in the Galactic disk, particularly for the isotopes 6Li, 9Be, and 10B, are well explained as being products of spallation processes in which CNO atoms in the interstellar medium (ISM) are broken up into LiBeB by collisions with CR protons or $`\alpha `$ particles (Reeves, Fowler & Hoyle 1970, Meguzzi, Audouze & Reeves 1971, Walker, Mathews & Viola 1985). In the last decade, extensive observations of LiBeB in population II, metal poor halo stars (MPHS) have turned up new and unexpected results, spurring controversy as to what type of CR sources and production mechanisms were operating in the halo of the early, forming Galaxy (see review by Vangioni-Flam, Cassé & Audouze 2000).

To date, most models of light element evolution in the early Galaxy have focused on strong shocks driven by supernovae (SNe) as the principal sources of CRs. Although the assumed CR composition, energy spectrum and the manner of LiBeB production vary in different models, a general consensus is that Be and B in MPHS mainly originate from the “inverse” spallation process, whereby CR CNO particles are transformed in flight into LiBeB by impinging on ISM H or He atoms (Duncan, Lambert & Lemke 1992). This can be realized if a sizable fraction of the CRs responsible for spallation comprise fresh, CNO-rich SN ejecta (e.g. Cassé, Lehoucq & Vangioni-Flam 1995, Vangioni-Flam et al. 2000, Ramaty et al. 2000, Parizot & Drury 1999, Suzuki, Yoshii & Kajino 1999, hereafter SYK, Suzuki & Yoshii 2001, hereafter SY), as opposed to CRs injected from the average ISM (e.g. Fields & Olive 1999).

The origin of 6Li in MPHS, which has been detected in only 3 stars so far (Hobbs 2000, Nissen 2000 and references therein), is more mysterious, as current models involving SN CRs face some difficulties. A peculiar aspect of Li is that in addition to spallation, the fusion process of CR $`\alpha `$ particles with ambient He atoms can be effective, and should actually dominate Li production at low metallicities. (Note that while both 7Li and 6Li are synthesized in comparable amounts, the CR-produced 7Li component is generally overwhelmed by the “Spite plateau” from primordial nucleosynthesis in the metallicity range under consideration; e.g. Ryan et al. 2001.) If the CR energy spectrum is taken to be a standard power-law distribution in momentum (§3), one requires a CR injection efficiency much higher than normally inferred to reproduce the 6Li observations, whether the CR composition is metal-enriched or not (Ramaty et al. 2000; SY, see their fig.2). Faring better are some models which envision CR acceleration by multiple shocks inside superbubbles, where an additional, $`\alpha `$-enriched low energy CR component with a hard spectrum and an energy cutoff at few 100 MeV is present (e.g. Vangioni-Flam et al. 1999, Parizot & Drury 1999). However, such a CR component is at the moment only hypothetical, without any observational support. This raises the question of whether there may have been other sources for 6Li.

This work investigates a new and more natural 6Li production scenario based on a previously unconsidered CR source: CRs accelerated at structure formation shocks, i.e. gravitational virialization shocks driven by the infall and merging of sub-Galactic gas clumps during the hierarchical build-up of structure in the early Galaxy. Such shocks are inevitable consequences in the currently standard theory of hierarchical structure formation. We show below that this picture gives a better explanation of the present data, and also provides a number of testable predictions for future observations that are quite distinct from SN CR models. Our scenario also embodies unique and important implications for understanding the formation of our Galaxy.

## 2. Cosmic Ray Sources in the Early Galaxy

SNe are known to release $`E_{\mathrm{SN}}10^{51}`$ erg of kinetic energy in each explosion, driving strong shocks into the ambient medium. These SN shocks are favorable sites for efficient CR acceleration through the first order Fermi mechanism (e.g. Blandford & Eichler 1987). A plausible value for the injection efficiency $`\xi _{\mathrm{SN}}`$, i.e. the fraction of the SN kinetic energy imparted to CRs, is $`1020\%`$, deduced from comparison of the SN rate and the energy content of CRs currently observed in the Galactic disk (Berezinskii et al. 1990, Chap.I, §4).

The global energetics of SNe for the early Galactic halo is more uncertain, but may be estimated as follows. We assume that the rate of SNe in the early Galaxy (dominated by core collapse SNe) was proportional to the global rate inferred from the cosmic star formation rate (e.g. Madau, della Valle & Panagia 1998). This is $`11.5`$ orders of magnitude greater at redshift $`z2`$, approximately the main epoch of halo star formation, compared to $`z0`$; with a present SN rate $`1/(30100\mathrm{y}\mathrm{r})`$, the average rate in the early halo is roughly $`1/(3\mathrm{y}\mathrm{r})`$. Taking the duration of active star formation in the halo to be $`5\times 10^8`$ yr (§3), the total number of SNe during the halo epoch is $`N_{\mathrm{SN}}1.5\times 10^8`$, and its integrated energy output is $`_{\mathrm{SN}}N_{\mathrm{SN}}E_{\mathrm{SN}}1.5\times 10^{59}`$ erg. If this is spread over gas with total mass $`M_g3\times 10^{11}\mathrm{M}_{}`$, the average specific energy input from SNe is $`ϵ_{\mathrm{SN}}_{\mathrm{SN}}\mu m_p/M_g0.15`$ keV per particle, where $`\mu 0.6`$ is the mean molecular weight. Alternatively, $`ϵ_{\mathrm{SN}}`$ can be gauged without recourse to $`M_g`$ from the total amount of heavy elements ejected by halo SNe. Taking for example oxygen, its abundance at the end of the halo phase should correspond to \[O/H\] $`1`$ or $`0.1\%`$ of the total gas mass $`M_g`$, implying a total mass in oxygen of $`10^4M_g`$. According to SN nucleosynthesis models (e.g. Nomoto et al. 1997), the oxygen yield per SN averaged over the mass distribution of progenitor stars is $`2M_{}`$. This leads to $`N_{\mathrm{SN}}5\times 10^5M_g`$ and $`ϵ_{\mathrm{SN}}=N_{\mathrm{SN}}\times 10^{51}\mathrm{erg}/M_g0.15`$ keV per particle, which is independent of $`M_g`$ and completely consistent with the above estimate.

SNe are not the only sources of mechanical energy (and hence CRs through shock acceleration) that may have been active in the early Galaxy. In the framework of the currently successful picture of hierarchical structure formation in the universe, large scale objects such as galaxies and clusters are formed through the merging and virialization of smaller subsystems, driven by gravitational forces acting on the dark matter. For each merging hierarchy, shocks should inevitably arise in the associated baryonic gas component, whereby the kinetic energy of infall is dissipated and the gas heated to the virial temperature of the merged halo (e.g. White & Rees 1978, Blumenthal et al. 1984, Peacock 1999). It is quite plausible that such structure formation (SF) shocks also accelerate CRs (e.g. Miniati et al. 2001).

A guide to how structure formation may have proceeded in our Galaxy, particularly for the Galactic halo, may be offered by the recent numerical simulations of Galaxy formation by Bekki & Chiba (2000, 2001, hereafter BC; see also Steinmetz & Müller 1995). The initial state is an overdense sphere in the expanding universe, superimposed on which are small scale fluctuations with a cold dark matter (CDM) power spectrum. Through gravitational instability, the perturbations on the smallest scales first develop into nonlinear sub-Galactic clumps, whose typical masses are expected to be $`10^610^7\mathrm{M}_{}`$. The clumps thereafter grow progressively more massive by merging with each other, until by a redshift of a few, most of the mass reside in the two largest clumps. These then merge into a single entity in the central region at redshift $`z2`$ (corresponding to the epoch of ‘formation’ or ‘final major merger’; e.g. Lacey & Cole 1993), whereby the majority of the infall kinetic energy is virialized. Stars eventually constituting the metal-poor halo are formed inside and dispersed from the merging clumps throughout the course of these events.

The energy dissipated at the main SF shock, that accompanying the major merger of the two most massive clumps, can be evaluated from the virial temperature of the merged system. The total mass of this system $`M_t`$ should be close to the total mass of the Galaxy $`M_{t,0}`$ measured today, as subsequent mass increase through accretion and mergers of external systems cannot have been very significant; otherwise, the thin disk, which must have formed slowly after halo formation and continues to exist today, would have been strongly disrupted (BC, Chiba & Beers 2001 and references therein). The virial temperature $`T_v`$ for a halo of total mass $`M_t`$ virializing at redshift $`z`$ is

$`k_BT_v={\displaystyle \frac{\mu m_pGM_t}{2r_v}}`$ $``$ $`0.26\mathrm{keV}\left({\displaystyle \frac{M_t}{3\times 10^{12}(h/0.7)^1\mathrm{M}_{}}}\right)^{2/3}`$ (1)
$`\times {\displaystyle \frac{1+z}{3}}{\displaystyle \frac{f(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },z)}{f(0.3,0.7,2)}}`$
where $`r_v`$ is the virial radius, and $`f(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },z)`$ is a factor which depends on the cosmology (given in e.g. Barkana & Loeb 2001). The numerical value corresponds to a lambda cosmology with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`h=0.7`$. The specific energy dissipated at the main SF shock is thus $`ϵ_{\mathrm{SF}}3k_BT_v/20.4`$ keV per particle, higher than the above estimate for SNe by a factor of $`2.6`$, if $`M_t=3\times 10^{12}\mathrm{M}_{}`$ is adopted. The CR contribution from SF shocks compared to SNe should similarly be higher, as the SF shock injection efficiency $`\xi _{\mathrm{SF}}`$ should not be too different from that for SNe (Miniati et al. 2001). We note that for $`M_{t,0}`$ ($`M_t`$), Wilkinson & Evans (1999) give $`M_{t,0}=1.9_{1.7}^{+3.6}\times 10^{12}M_{}`$ utilizing kinematic information of satellite galaxies and globular clusters, while Zaritsky (1999) puts a lower bound of $`M_t1\times 10^{12}M_{}`$ by summarizing a variety of methods. The true value probably lies inside $`M_t(15.5)\times 10^{12}M_{}`$, corresponding to a range $`ϵ_{\mathrm{SF}}/ϵ_{\mathrm{SN}}1.34.0`$. (The recent analysis of Sakamoto, Chiba & Beers (2002) using a large sample of mass tracers results in $`M_t=2.5_{1.0}^{+0.5}\times 10^{12}M_{}`$ if the satellite Leo I is bound to the Milky Way, and $`M_t=1.8_{0.7}^{+0.4}\times 10^{12}M_{}`$ if not.) CRs accelerated by SF shocks should therefore be at least as important as SN CRs, and may well dominate at early epochs.

Another, more speculative possibility is outflows powered by massive black hole(s), which may have been active in the early Galaxy (e.g. Silk & Rees 1998). Our Galactic center today is known to harbor a black hole of mass $`M_{\mathrm{BH}}3\times 10^6\mathrm{M}_{}`$ (Melia & Falcke 2001), and if a fraction $`\eta `$ of its rest mass energy could have been converted into kinetic energy of a jet or a wind, the total ejected energy would amount to $`_{\mathrm{BH}}\eta M_{\mathrm{BH}}c^25.4\times 10^{59}\mathrm{erg}(\eta /0.1)(M_{\mathrm{BH}}/3\times 10^6\mathrm{M}_{})`$. Strong shocks induced by this outflow may then provide CRs energetically comparable to or greater than SNe or SF (c.f. Crosas & Weisheit 1996, SY). Black hole outflows are thus potentially interesting; however there are numerous ambiguities with such a picture, and this paper will concentrate on CRs from SF and SNe.

## 3. Model Description

SYK and SY developed a model of light element evolution based on SN CRs, taking into account the inhomogeneous nature of SN-induced chemical evolution in the early Galactic halo (Tsujimoto, Shigeyama & Yoshii 1999, hereafter TSY), and including time-dependent calculations of CR propagation. This model is extended in order to incorporate the effect of CRs from SF shocks. Chemical evolution of the halo initiates with the first generation stars forming at time $`t=0`$, and lasts until $`t=0.6`$ Gyr when SN remnants have swept up all the gas in the halo and \[Fe/H\] has reached $`1.5`$.

We employ assumptions and parameters deemed most plausible for the SN CRs. The CR injection efficiency is taken to be $`\xi _{\mathrm{SN}}=0.15`$ for each SN of kinetic energy $`E_{\mathrm{SN}}=10^{51}\mathrm{erg}`$. A single power-law distribution in particle momentum is adopted for the CR spectrum, the form expected from standard shock acceleration theory (Blandford & Eichler 1987); in terms of particle energy per nucleon $`E`$, this is proportional to $`(E+E_0)[E(E+2E_0)]^{\frac{\gamma _{\mathrm{SN}}+1}{2}}`$, where $`E_0=930`$MeV is the nucleon rest mass energy. The injection spectral index is chosen to be $`\gamma _{\mathrm{SN}}=2.1`$, appropriate for strong SN shocks, and consistent with the source spectrum inferred for present-day CRs. As assumed in SYK and SY, the composition of SN CRs is a mixture of SN ejecta containing freshly synthesized CNO and Fe and the ambient ISM swept up by the SN blastwave. This is quantified by the parameter $`f_{\mathrm{CR}}`$, the fractional mass of the ISM in the SNR shell relative to the SN ejecta mass that goes into CR acceleration. Since SNe with different progenitor masses give varying nucleosynthetic yields (e.g. TSY), the composition of CRs propagating in the ISM also corresponds to an ensemble of SN yields. The total injected SN CR flux at time $`t`$ is normalized by the SN rate at that epoch, which is specified by our chemical evolution model, along with the above values for $`\xi _{\mathrm{SN}}`$ and $`E_{\mathrm{SN}}`$. In order to reproduce the observed abundances of Be and B with our calculated results, the parameter $`f_{\mathrm{CR}}`$ is uniquely determined to be $`2\times 10^4`$ (see fig.2 in SY). This then fixes the ratio of $`\alpha `$’s to protons in SN CRs to be $`0.14`$ (the $`\alpha `$/p ratio is $`0.08`$ for the ISM and $`0.2`$ for the SN ejecta, SY), which in turn fixes the resultant 6Li yield.

CRs from SF shocks are markedly different from SN CRs in a number of important ways. First and foremost, SF shocks do not synthesize fresh CNO nor Fe, so that the composition of these CRs is completely ascertained by the pre-existing ISM. When the ISM is metal-poor, these shocks induce very little Be or B production through inverse spallation, and are only efficient at spawning Li via $`\alpha \alpha `$ fusion. Second, SF shocks are not necessarily strong ones, particularly for major mergers of systems with comparable masses (Miniati et al. 2001), which applies to the case mainly envisaged here (§2). Shock acceleration should then lead to injection indices $`\gamma _{\mathrm{SF}}`$ steeper than the strong shock limit value of 2, which works in favor of Li production (§4). Major merger shocks may possess Mach numbers as low as $`23`$, corresponding to $`\gamma _{\mathrm{SF}}2.53.3`$; $`\gamma _{\mathrm{SF}}=3`$ is generally chosen below. As with SN CRs, we take the spectral shape to be a momentum power-law distribution and the injection efficiency to be $`\xi _{\mathrm{SF}}=0.15`$.

A further distinction from SNe is that the SF CR flux should not entail any direct dependence on the metallicity, which is in fact an obstacle to predictive modeling. While the hierarchical growth of structure with respect to redshift or cosmic time may be evaluated in concrete ways using certain formalisms (e.g. Lacey & Cole 1993), relating this to \[Fe/H\] requires additional knowledge of how the combination of star formation, SN nucleosynthesis and chemical evolution in our Galaxy proceeded with redshift or cosmic time. This involves large uncertainties, and is not specified in our chemical evolution model. As a first step, we choose to describe the time evolution of SF CRs in a simple, parameterized way, assuming a ‘step function’ behavior: SF CRs begin to be injected from a certain time $`t_{\mathrm{SF}}`$, maintains a constant flux for a duration $`\tau _{\mathrm{SF}}`$, and then returns to zero. We take the injection duration $`\tau _{\mathrm{SF}}`$ to be roughly the dynamical time of the major merger, $`3\times 10^8`$ yr. The injection is also assumed to be uniform, since the effect of the main SF shock should be global throughout the gas under consideration, unlike SNe. The injected SF CR flux integrated over $`\tau _{\mathrm{SF}}`$ is normalized to the above values for $`\xi _{\mathrm{SF}}`$ and $`ϵ_{\mathrm{SF}}`$. The true evolutionary behavior of SF shock activity relative to metallicity should actually be probed through future observations of 6Li at low \[Fe/H\] (§4, §5).

After injection by either SN or SF shocks, the spectral flux $`F_i(E,t)`$ for each CR element $`i`$ evolves with time during subsequent interstellar propagation. This is obtained from time-dependent solutions of the CR transport equation for a leaky box propagation model (Meneguzzi et al. 1971, SYK, SY), including escape and losses due to nuclear destruction and ionization. At $`E`$ 50 - 500 MeV per nucleon, the most relevant energies for LiBeB synthesis, the effective loss length due to nuclear destruction $`\mathrm{\Lambda }_{\mathrm{n},i}20\mathrm{g}\mathrm{c}\mathrm{m}^2`$ for CR $`\alpha `$’s and CNOs (Malaney & Butler 1993), and the ionization loss length $`\mathrm{\Lambda }_{\mathrm{ion}}(150)(A/Z^2)\mathrm{gcm}^3`$ for CR particles with mass number $`A`$ and charge $`Z`$ (Northcliffe & Schilling 1970).

The escape length in the early halo is very uncertain, but we may speculate from the presumed properties of the gas in the merged system following the main SF shock, which could be of total mass $`M_g3\times 10^{11}\mathrm{M}_{}`$ and spatial extent $`R`$ 10 kpc (BC). If the CR diffusion coefficient in this medium is assumed to be similar to the value inferred for the halo today, $`D_{\mathrm{CR}}10^{29}\mathrm{cm}^2\mathrm{s}^1`$ (e.g. Berezinskii et al. 1990, chap.3, §3), the escape length is roughly $`\mathrm{\Lambda }_{\mathrm{esc}}(3M_g/4\pi R^3)(R^2/6D_{\mathrm{CR}})c200\mathrm{g}\mathrm{c}\mathrm{m}^2`$. To be compared are $`\mathrm{\Lambda }_{\mathrm{n},i}`$ and $`\mathrm{\Lambda }_{\mathrm{ion}}`$ above, implying that escape is unimportant here, and ionization and nuclear destruction losses determine the CR spectra in the energy range effective for light element production. We set $`\mathrm{\Lambda }_{\mathrm{esc}}=100\mathrm{g}\mathrm{c}\mathrm{m}^2`$ in our calculations below, although the exact value is irrelevant so long as this is larger than $`\mathrm{\Lambda }_{\mathrm{ion}}`$ or $`\mathrm{\Lambda }_{\mathrm{n},i}`$. It is emphasized that under such conditions, the LiBeB production rate does not depend on the gas density $`n_g`$ (as opposed to the situation considered in e.g. Fields et al. 2001): since both ionization and nuclear destruction loss rates are proportional to $`n_g`$, any increase (decrease) in the density of gas containing target particles is always compensated by a decrease (increase) in the propagated CR flux due to increased (decreased) losses.

Using the transported spectra, we calculate the CR production of LiBeB in the ISM including all three types of reactions described in $`\mathrm{\S }1`$: forward spallation of ISM CNO by CR protons and $`\alpha `$’s, inverse spallation of CR CNOs by ISM H and He, and the fusion of CR $`\alpha `$’s with ISM He. For more details, consult SYK and SY. We do not consider here the potentially complicating effects of stellar depletion, which are highly uncertain at the moment (e.g. Pinsonneault et al. 1999) and await further studies.

## 4. Results and Discussion

For our calculations, we have selected the following sets of parameters for $`t_{\mathrm{SF}}`$, $`\tau _{\mathrm{SF}}`$ and $`\gamma _{\mathrm{SF}}`$, respectively, labeled models I - V: I (0.12, 0.1, 3), II (0.22, 0.1, 3), III (0.32, 0.1, 3), IV (0.22, 0.1, 2) and V (0.1, 0.5, 3), where $`t_{\mathrm{SF}}`$ and $`\tau _{\mathrm{SF}}`$ are in units of Gyr; model VI is the case of SN CRs only. These were chosen to provide results exemplary of light element production by SF CRs, in contradistinction to that by SN CRs. The evolution of 6Li and Be vs. metallicity calculated for each model until the end of halo chemical evolution (\[Fe/H\]$`1.5`$) is shown in Fig.1, along with the current observational data for 6Li (Smith et al. 1998, Cayrel et al. 1999, Nissen et al. 2000) and Be (Boesgaard et al. 1999) in MPHS.

We discuss some salient points regarding these results. First, it is confirmed that with our fiducial parameters, production by SN CRs alone (VI) works very well for the observed Be (and B, not shown), yet falls short of the observed 6Li. Accounting for this by SN CRs demands a much larger value of $`E_{\mathrm{SN}}`$, $`\xi _{\mathrm{SN}}`$ and/or $`\gamma _{\mathrm{SN}}`$, which must be supplemented with a larger $`f_{\mathrm{CR}}`$ in order for Be and B to be consistently reproduced together (see fig.2 in SY).

In contrast, with reasonable values for $`ϵ_{\mathrm{SF}}`$, $`\xi _{\mathrm{SF}}`$ and $`\gamma _{\mathrm{SF}}`$, production by SF CRs is capable of explaining the current 6Li data quite adequately. This mainly owes to two facts: 1) SF shocks are more energetic than (or at least as energetic as) SN shocks, as estimated in §2, and 2) SF CRs can generate 6Li at early epochs independently of the metallicity. Regardless of the early evolutionary behavior, identical 6Li abundances are attained at the end of the halo phase for a given $`\gamma _{\mathrm{SF}}`$ (I - III, V), since this is determined by the time-integrated CR flux, for which we had assumed a fixed value. Compared to a flat spectral index of $`\gamma _{\mathrm{SF}}=2`$ (IV), a steeper one of $`\gamma _{\mathrm{SF}}=3`$, more appropriate for low Mach number SF shocks (§3), results in a larger 6Li yield, by about factor of 3. This is because with a constant total CR energy, a steeper index implies a larger CR flux in the subrelativistic energy range $`E10`$$`100`$ MeV, where the 6Li production cross section peaks (see fig.1 in SY). Note that a conservatively lower specific energy for SF shocks, $`ϵ_{\mathrm{SF}}0.15`$ keV/particle (i.e. comparable to SNe), can still be consistent with the available 6Li data provided that $`\gamma _{\mathrm{SF}}3`$.

As already stressed, 2) is a consequence of 6Li synthesis being dominated by $`\alpha \alpha `$ fusion, and SF shocks not creating any new Fe. Depending on the onset time and duration of the SF shock, the 6Li abundance may potentially reach large values quickly at very low metallicity, which can be followed by a plateau or a very slow rise. On the other hand, SNe unavoidably give forth to freshly synthesized Fe, so a correlation with Fe/H must arise; in fact 6Li/H vs. Fe/H for SN CRs can never be much flatter than linear (e.g. Vangioni-Flam et al. 1999, Parizot & Drury 1999, Fields & Olive 1999, Ramaty et al. 2000). Moreover, since SF shocks do not eject fresh CNO, they produce very little Be or B through spallation; only a minuscule contribution can appear as the ISM becomes metal-enriched. This may allow an extremely large 6Li/Be ratio at low \[Fe/H\]. Conversely, we see that SN CRs must play an indispensable role in generating the Be and B observed in MPHS.

Independent of the particular evolutionary parameters, the following abundance trends are characteristic of SF CR production and should serve as distinguishing properties of the scenario for future observations. Going from high to low metallicity: a plateau or a very slow decrease in $`\mathrm{log}^6`$Li/H vs. \[Fe/H\], followed by a steeper decline in some range of \[Fe/H\] corresponding to the main epoch of SF; a steady increase in 6Li/Be, possibly up to values exceeding $`100`$, also followed by a downturn. These traits are very distinctive and not expected in SN CR models, for which the slope of $`\mathrm{log}^6`$Li/H - \[Fe/H\] must be $`1`$ or greater, and the 6Li/Be ratio constant at sufficiently low \[Fe/H\]. Distinction from any production processes in the early universe (e.g. Jedamzik 2000) should also be straightforward, as they predict a true plateau down to the lowest \[Fe/H\], in contrast to an eventual decrease for SF CR models. Further, unique diagnostic features are discussed in §5.

Observing 6Li in MPHS has proven to be difficult in the past, as measurement of its weak isotopic shift feature relative to the much stronger 7Li line requires spectroscopy with very high resolution and signal to noise. The present database of only 3 positive detections in a narrow range of \[Fe/H\] ($``$ -2.5 – -2) plus a number of upper limits (Hobbs 2000, Nissen 2000 and references therein) is obviously insufficient for distinction between our SF CR model and various other models. However, the new generation of large aperture telescopes equipped with high resolution spectrographs, e.g. the Subaru HDS or VLT/UVES, should in the near future bring about a larger sample of MPHS with accurately determined 6Li/H and 6Li/Be over a wide range of \[Fe/H\], and greatly help us toward deciphering the true origin of 6Li in MPHS.

## 5. 6Li as Fossil Record of Dissipative Processes during Galaxy Formation

A truly unique and intriguing aspect of the SF CR picture is that 6Li in MPHS can be interpreted and utilized as a fossil record of dissipative gas dynamical processes in the early Galaxy. Of particular interest are various correlations expected between the 6Li abundance and the kinematic properties of the stars. On the one hand, 6Li arises as a consequence of gaseous dissipation through gravitational shocks, and survives to this day as signatures of the dynamical history of hierarchical structure formation in the early Galaxy. On the other, the kinematic characteristics of stars presently observed should reflect the past dynamical state of their parent gas systems, because once stars form, they become collisionless and have long timescales for phase space mixing (e.g. Chiba & Beers 2000). Interesting relationships may then exist among the two observables.

For example, recent intensive studies of the structure and kinematics of MPHS in our Galaxy based on Hipparcos data (e.g. Chiba & Beers 2000, 2001) have elucidated the detailed characteristics of our Galaxy’s halo, such as its two-component nature: an inner halo which is flattened and rotating, and an outer halo which is spherical and non-rotating. This dichotomy has been suggested to result from differences in the physical processes responsible for their formation, dissipative gas dynamics being crucial for the former, and dissipationless stellar dynamics determining the latter. The numerical simulations of Galaxy formation by BC support this conjecture (see §2 for a rough outline): the outer halo forms through dissipationless merging of small sub-Galactic clumps that have already turned into stars (c.f. Searle & Zinn, 1978), whereas the inner halo mainly forms through dissipative merging and accretion of larger clumps that are still gas rich (c.f. Eggen, Lynden-Bell & Sandage, 1962). In our scenario, 6Li production is a direct outcome of the principal gas dissipation mechanism of gravitational shock heating. If the above inferences on the formation of halo structure are correct, 6Li should be systematically more abundant in stars belonging to the inner halo compared to those of the outer halo, a clearly testable prediction. An important prospect is that 6Li may provide a quantitative measure of the effectiveness of gas dynamical processes during formation of halo structure, rather than just the qualitative deductions allowed by kinematic studies.

Another possibility regards the main epoch of SF with respect to metallicity (§3). As already mentioned, the relation between 6Li/H and Fe/H should mirror the time evolution of dissipative energy release through SF shocks, but is complicated by being convolved with the uncertain ingredients of star formation and chemical evolution. The Fe abundance can be a bad tracer of time, especially at low \[Fe/H\] where the effects of dispersion in SN yields can be extremely large (TSY, SY). Stellar kinematics information may offer a handle on this problem. The observed relation between \[Fe/H\] and $`<V_\varphi >`$, the mean azimuthal rotation velocity of MPHS, seems to manifest a distinctive kink around \[Fe/H\] $``$ -2 (Chiba & Beers 2000). Through chemo-dynamical modeling of the early Galaxy, BC have proposed that this kink may correspond to the epoch of the major merger (§3). If this was true, a simple expectation in the context of the SF CR model is that 6Li/H should be just rising near this value of \[Fe/H\], which is the range occupied by the currently 6Li-detected stars; also expected are a steep decline at lower \[Fe/H\], as well as a plateau or slow rise at higher \[Fe/H\] (i.e. close to model III in fig.1). However, any inferences related to \[Fe/H\] are always subject to the chemical evolution ambiguities. A more reliable and quantitative answer may be achieved by looking for correlations between 6Li/H and $`<V_\varphi >`$ without recourse to \[Fe/H\], as 6Li is a direct and pure indicator of dynamical evolution in the early Galaxy.

Thus the SF CR model for 6Li bears important implications for understanding how our Galaxy formed. If the above mentioned trends are indeed observed, it would not only confirm the structure formation origin of 6Li, but may potentially point to new studies of “6Li Galactic archaeology”, whereby extensive observations of 6Li in MPHS can be exploited as a robust and clear-cut probe of dissipative dynamical processes that were essential for the formation of the Galaxy.

## 6. Summary

We have put forth a new scenario for the currently puzzling origin of 6Li observed in metal-poor halo stars: production by cosmic rays accelerated at structure formation shocks, driven by the hierarchical infall and merging of sub-Galactic structure during the formation of our Galaxy. Several predictions are quite distinct from models involving supernova cosmic rays and clearly testable in the near future, such as the behavior of 6Li/H vs. \[Fe/H\] and 6Li/Be vs. \[Fe/H\] at low \[Fe/H\], as well as possible correlations between 6Li/H and the kinematic properties of halo stars. Since 6Li can be construed as a fossil record of dissipative processes during Galaxy formation, further observations of this isotope in halo stars may offer us a unique and invaluable insight into the past dynamical history of our Galaxy.

The authors are very grateful to M. Chiba for extensive and fruitful discussions. We also thank G. Mathews, T. Kajino and M. Nagashima for valuable conversations, Y. Yoshii for assistance with our code, and the anonymous referee for helpful comments. T.K.S. is supported by the JSPS Research Fellowship for Young Scientists, grant 5936.