text:
# Ly𝛼 Transfer in a thick, dusty, and static medium

## 1 Introduction

Recent advances of large telescopes and CCD techniques have been unveiling very remote astrophysical objects such as primeval star forming galaxies and damped Ly$`\alpha `$ galaxies. The Ly$`\alpha `$ emission profiles of these primeval galaxies are categorized into three types : (1) pure Ly$`\alpha `$ emission, (2) asymmetric or P Cyg type one, (3) broad absorption in damping wings. Moreover, local starburst galaxies also show similar characteristics (e.g. Kunth et al. 1998). We will not list up the full inventory of the starburst galaxies mentioned above, but the issue was discussed in our previous paper (Lee & Ahn 1999). It is plausible that dusts play an important role in the formation of the Ly$`\alpha `$ line profile. Since Ly$`\alpha `$ sources in a starburst galaxies are usually embeded in optically thick and dusty media, we must investigate the effects of dusts on the radiative transfer and line formation of Ly$`\alpha `$ in the media.

Recently Meurer et al. (1999) studied the dust extinction in local starburst galaxies, and applied the result to the primeval galaxies in the Hubble Deep Field. They found a tight relationship between the $`\beta `$ indexes of ultraviolet (UV) continua and the ratios of far-infrared (FIR) to UV fluxes of the galaxies. The relationship is established because dusts absorb Ly$`\alpha `$ line photons and UV continua of $`\lambda _0>912\AA `$, and then re-emit in FIR. Leitherer et al. (1999) suggested that the Ly$`\alpha `$ luminosity is as large as $`10\%`$ of the total UV luminosity in the star forming regions. Hence it is necessary to know what fraction of Ly$`\alpha `$ photons is subject to dust extinction during the transfer in a very thick neutral medium.

Studies of the Ly$`\alpha `$ line transfer in an optically thick and static medium has a long history. Unno (1955) formulated the Ly$`\alpha `$ line transfer in a dust-free medium, and Osterbrock (1962) proposed a simple physical picture for understanding the resonance line transfer in a thick medium. Adams (1972) revised Osterbrock’s picture by a Monte Carlo method and gave a heuristic explanation on the problem. In an analytical way, Harrington (1973) solved the problem, for which Neufeld (1990) also gave a more general solution. However, the studies thus far are limited to the cases where diffusion approximation is valid and little attention has been paid to the polarization of the Ly$`\alpha `$ flux emergent from those media with anisotropic geometry and/or kinematics. More realistic treatment of the Ly$`\alpha `$ line transfer is clearly necessary.

We have developed a sophisticated Monte Carlo code that describes the radiative trasfer of Ly$`\alpha `$ resonnance photons in optically thick, dusty, and static media, ranging from moderate optical depths to extremely large optical depths. In Section II, we describe the basic theory including the numerical method of our code. The basic results are presented in Section III, and the possible applications of our code are described in Section IV.

## 2 Basic Theory

### 2.1 Configuration and Review

The medium used in our investigation is plane-parallel, in which hydrogen atoms and dust particles are mixed homogeneously. In Fig. 1, we show the configuration. We consider the case that a source is located in the mid-plane and radiates photons isotropically. However, it is easy to extend to other configurations.

Osterbrock (1962) proposed a simple physical picture based on a simplified treatment of random walk processes, and predicted that $`N\tau _0^2`$, where $`N`$ is the mean number of scattering and $`\tau _0`$ is the line center optical depth. But Adams (1972) solved the problem by a Monte Carlo method and showed that $`N\tau _0`$ rather than $`N\tau _0^2`$, where $`a`$ is the Voigt parameter. Also $`N`$ and the peak of emergent spectra is dependent upon $`a\tau _0`$ rather than $`\tau _0`$ alone. He made an explanation that the line transfer process is dominated by wing scattering, and is well approximated by random walk processes. He also pointed out that the transfer of Ly$`\alpha `$ photons in a thick medium can be approximated as a diffusion process both in real space and in frequency space. He suggested that in moderately thick media (i.e. $`10^3<\tau _0<10^3/a`$) Ly$`\alpha `$ photons escape by a single longest flight, and for extremely thick media (i.e. $`\tau _0>10^3/a`$) by a single longest excursion.

Harrington (1973) solved the radiative transfer equation presented by Unno (1955) in an analytic manner. He adopted the diffusion approximation and used the Sturm-Liouville theory to obtain the analytic solution. However, he considered only for one specific configuration, in which a monochromatic ($`\lambda _0=1216\AA `$) source is located at the center of a plane-parallel medium. Neufeld (1990) gave analytic solutions for a more generalized version of the problem, where the source produces arbitrary initial frequencies and various dust extinctions were included.

In the present paper, we will consider the same problems with more general optical depths using Monte Carlo techniques, and the results obtained by Neufeld (1990) will be used to assure the validity of our method. However, it is noted that our Monte Carlo approach is more faithful to the problem by resorting to less approximation than other approaches mentioned above.

### 2.2 Optical Depth

In the configuration stated in the previous subsection, the optical depth is given by

$`\tau _\nu (s)={\displaystyle _0^s}𝑑l{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑v_zn(v_z)\sigma _\nu .`$ (1)
In a static medium, the Ly$`\alpha `$ scattering cross section is given by

$`\sigma _\nu ={\displaystyle \frac{\pi e^2}{m_ec}}f_{12}{\displaystyle \frac{\mathrm{\Gamma }/4\pi ^2}{(\nu \nu _0)^2+(\mathrm{\Gamma }/4\pi )^2}},`$ (2)
where $`m_e`$ is the electron mass, $`c`$ the light velocity, $`f_{12}=0.4126`$ the oscillator strength for hydrogen Ly$`\alpha `$, $`\mathrm{\Gamma }=A_{21}=3.1\times 10^9\mathrm{s}^1`$ the damping constant, and $`\nu _0`$ the frequency of the line center. When we include the thermal motions of scatterers, $`\nu `$ is tranformed by $`\nu ^{}=\nu \nu _0\frac{v_z}{c}`$, and the thermal motion is described by the Maxwellian distribution,

$`n(v_z)dv_z=n_{HI}{\displaystyle \frac{1}{\sqrt{\pi }v_{th}}}e^{(\frac{v_z}{v_{th}})^2}dv_z.`$ (3)
Here $`v_{th}^2=2k_BT_g/m_H`$, where $`m_H`$ is the hydrogen mass, $`k_B`$ the Boltzmann constant, $`n_H`$ the neutral hydrogen density, and $`T_g`$ the gas temperature.

Substituting Eqs.(2) and (3) into Eq.(1), we obtain the expression for the optical depth,

$`\tau _x(s)`$ $`=`$ $`{\displaystyle \frac{\pi ^{1/2}e^2}{m_ec}}f_{12}n_{HI}s{\displaystyle \frac{1}{\nu _D}}{\displaystyle \frac{a}{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑u{\displaystyle \frac{e^{u^2}}{(xu)^2+a^2}}`$ (4)
$`=`$ $`1.41\times 10^{13}T_{g4}^{1/2}n_{HI}sH(x,a),`$ (5)
where $`T_{g4}T_g/10^4\mathrm{K}`$,

$`\nu _D={\displaystyle \frac{v_{th}}{c}}\nu _0,`$ (6)
$`x{\displaystyle \frac{\mathrm{\Delta }\nu }{\mathrm{\Delta }\nu _D}}{\displaystyle \frac{\nu \nu _0}{\mathrm{\Delta }\nu _D}},`$ (7)
$`u{\displaystyle \frac{v_z\nu _0}{c\mathrm{\Delta }\nu _d}},`$ (8)
and

$`a{\displaystyle \frac{\mathrm{\Gamma }}{4\pi \mathrm{\Delta }\nu _d}}.`$ (9)
Here $`H(x,a)`$ is the Hjerting function or the H-function, and the specific values in Rybicki & Lightman (1979) are used. Also the Voigt parameter $`a`$ is given by

$`a=4.7\times 10^4T_{g4}^{1/2}.`$ (10)
Therefore, the total optical depth for a given system is

$`\tau _x=1.41\times 10^{13}T_{g4}^{1/2}N_{HI}H(x,a)=\tau _0H(x,a),`$ (11)
where the column density of neutral hydrogen, $`N_{HI}=n_{HI}L`$ with $`L`$ being the physical size of the slab, and $`\tau _01.41\times 10^{13}T_{g4}^{1/2}N_{HI}`$ is the line center optical depth.

### 2.3 Monte Carlo Code

In this subsection, we present a detailed description of the Monte Carlo procedure. There are a few Monte Carlo approaches to the resonance line transfer in an optically thick and static medium in the literature (e.g. Adams 1972, Harrington 1973, Gould & Weinberg 1996).

The Monte Carlo code begins with the choice of the frequency and the propagation direction $`\stackrel{}{k}_i`$ of an incident photon from an assumed Ly$`\alpha `$ profile, which is assumed to be monochromatic in our study.

Then we determine the next scattering site separated from the initial point by the propagation length defined by

$`S\tau _0s={\displaystyle \frac{\tau }{H(x,a)+\tau _d/\tau _0}}{\displaystyle \frac{\tau }{H(x,a)}},`$ (12)
where the optical depth is assumed to be composed of hydrogen ($`\tau _0H(x,a)`$) and dust ($`\tau _d`$) parts. Here, $`\tau _d`$ is the dust optical depth which is treated as a constant for the emission line and whose Galactic value is given by Draine & Lee (1984), and $`\tau =\mathrm{ln}(𝐑)`$, where $`𝐑`$ is a uniform random number in the interval $`(0,1)`$ generated by a subroutine $`ran2()`$ suggested by Press et al.(1989).

The emitted photon traverses a distance $`S`$, and is scattered off by hydrogen atoms until $`|z|>\tau _0`$ for a slab geometry. In this scattering event the frequencies of the absorbed photon and the re-emitted one in the rest frame of the scatterer should be matched.

In contrast to the case of a thin medium, for a very thick medium the scattering in the damping wings is not negligible. A careful treatment needs to be exercised to distinguish the scattering in the damping wings from the resonance scattering, because they show quite different behaviors in the properties including the scattering phase function and the polarization (Lee & Blandford 1997).

Because the natural line width is much smaller than the Doppler width, the local velocity of the scatterers that can resonantly scatter the incident photon is practically a single value. However, when the scattering occurs in the damping wings, the local velocity of the scatterer may run a rather large range. Therefore, in order to enhance the efficiency of the Monte Carlo method, it is desirable to determine the scattering type before we determine the local velocity $`\stackrel{}{u}`$ of the scatterer.

We present a more quantitative argument about the preceding remarks. Under the condition that a given photon is scattered by an atom located at a position $`s`$, the local velocity component $`u`$ along the direction $`\stackrel{}{k}_i`$ is chosen from the normalized distribution,

$`f(u)={\displaystyle \frac{e^{u^2}}{(xu)^2+a^2}}\left[{\displaystyle \frac{\pi }{a}}H(x,a)\right]^1.`$ (13)
Here, the Hjerting function or the Voigt function $`H(x,a)`$ is evaluated by a series expansion in $`a`$, i.e.,

$`H(x,a)`$ $`=`$ $`H_0(x,a)+aH_1(x,a)+a^2H_2(x,a)`$ (14)
$`+`$ $`a^3H_3(x,a)+\mathrm{},`$
where $`H_n(x,a),n=0,1,2,3`$ are tabulated by Gray (1992).

Because of the smallness of $`a`$, the function $`f`$ has a sharp peak around $`ux`$, for which the scattering is resonant. Therefore, the probability $`P_r`$ that a given scattering is resonant is approximately given by

$`P_r`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}d(\mathrm{\Delta }u){\displaystyle \frac{e^{x^2}}{(\mathrm{\Delta }u)^2+a^2}}\left[{\displaystyle \frac{\pi }{a}}H(x,a)\right]^1`$ (15)
$`=`$ $`{\displaystyle \frac{e^{x^2}}{H(x,a)}}.`$
The probability that scattering occurs in the damping wings is

$`P_{nr}=1P_r.`$ (16)
In the code we determine the scattering type in accordance with the scattering type probabilities $`P_r`$ and $`P_{nr}`$. If scattering is chosen to be resonant, then we set $`u=x`$. Otherwise, the scattering occurs in the damping wings, and $`u`$ is chosen in accordance with the velocity probability distribution given by Eq.(12).

We give the propagation direction $`\stackrel{}{k}_f`$ of a scattered photon in accordance with the phase function completely faithful to the atomic physics (Lee & Ahn 1998). The scattered velocity component $`v_{}`$ perpendicular to the initial direction $`\stackrel{}{k}_i`$ on the plane spanned by $`\stackrel{}{k}_i`$ and $`\stackrel{}{k}_f`$ is also governed by the Maxwell-Boltzmann velocity distribution, which is numerically obtained using the subroutine $`\mathrm{𝑔𝑎𝑠𝑑𝑒𝑣}()`$ suggested by Press et al. (1989). The contribution $`\mathrm{\Delta }x`$ of the perpendicular velocity component $`v_{}`$ to the frequency shift along the direction of $`\stackrel{}{k}_f`$ is obviously

$`\mathrm{\Delta }x=v_{}[1(\stackrel{}{k}_i\stackrel{}{k}_f)^2]^{1/2}/c.`$ (17)
Therefore, the frequency shift $`x_f`$ of the scattered photon is given by

$`x_f=x_iu+u(\stackrel{}{k}_i\stackrel{}{k}_f)+v_{}[1(\stackrel{}{k}_i\stackrel{}{k}_f)^2]^{1/2},`$ (18)
where $`x_i`$ is the frequency shift of the incident photon.

In each scattering event the position of scattered photon is checked, and its path length is added. If the photon escapes from the medium or $`|z|<\tau _0`$, we collect that photon according to its frequency and escaping direction. In collecting photons, we add the weighted fraction considering the dust extinction corresponding to the total path length of the photon. The whole procedure is repeated until typically about $`10^3`$ photons in each frequency bin are collected.

## 3 Results

### 3.1 Tracing the Scattering Processes

Adams (1972) gave a physical picture that describes the Ly$`\alpha `$ line transfer in extremely thick media of $`\tau _0>10^3/a`$, where $`a=4.71\times 10^4`$ for $`T_g=10^4\mathrm{K}`$.

Initially a core photon becomes a wing photon after a sufficient number of scatterings and encountering with a violently moving atom. Due to the small optical depths, wing photons traverse much further in physical space than core photons. When $`\tau _0<10^3`$, a photon can directly escape from a medium once it becomes a wing photon. This process is called ‘a single longest flight’ by Adams (1972).

As the line center optical depth gets larger to be an intermediate optical depth, $`10^3<\tau _0<10^3/a`$, wing photons cannot directly escape from the medium, but experience a large number of core scatterings followed by a few wing scatterings. However, after these wing scatterings, photons may return back to the core because of the so called ‘the restoring force’ (Osterbrock 1962).

For the case of the extremely thick medium we can describe this process using a diffusion approximation, where the $`rms`$ frequency shift per scattering is the thermal Doppler width ($`v_{th}`$) or $`x=1`$, and the mean frequency shift per scattering is $`\mathrm{\Delta }x=1/|x|`$ that plays the role of the restoring force. This restoring force is caused by relatively large probability that photons scatter in the core. However, for a medium with an intermediate optical depth, the first variance of frequency is not quite large because the optical depth is not enough to enlarge the probability for wing scattering at frequencies far from the line center. Hence, photons often experience a small number of single longest flights before escape, and we may call these processes as ‘random wandering’. Therefore, in this particular case, we think that a simple diffusion approximation is not adequate and that a more accurate approach such as a Monte Carlo method is needed.

For a medium of an even higher optical depth, $`\tau _0>10^3/a`$, a series of wing scatterings occur when wing photons are generated. In this case the wandering occurs both in real space and in frequency space. This process is called ‘an excursion,’ and photons escape the medium by ‘a single longest excursion’ (Adams 1972).

These processes can be seen in our Monte Carlo results. In Fig. 2 we have shown the trajectory of a photon that transfers in a medium with $`\tau _0=10^4`$. At the 8055th scattering, the photon is scattered off by an atom that moves fast and became a wing photon. At the same time, the photon traverses a longer distance in physical space, and also wing scattering happens.

Before escaping from the medium, the photon experiences a series of wing scatterings. In Fig. 3 we show the history of the same photon, and here we can see that the last scattering occurs in the wing. We checked that the larger $`\tau _0`$, the larger the number of last wing scatterings, which may be regarded as ‘a single longest excursion.’ These last wing scatterings possibly induce polarization because wing scattering has a Rayleigh phase function (Stenflo 1980). We have conducted an investigation on this topic, and the results will be published in another paper.

### 3.2 Emergent Ly$`\alpha `$ Profile for Dust Free Media

Neufeld (1990) derived an analytic solution, for $`(a\tau _0)^{1/3}1`$. According to Unno (1955), Harrington (1973), and Neufeld (1990), for the cases of dust-free media, the Ly$`\alpha `$ line transfer equation is given by

$`{\displaystyle \frac{^2J}{\tau ^2}}+{\displaystyle \frac{^2J}{\sigma ^2}}=3\varphi (x){\displaystyle \frac{E(\tau ,x)}{4\pi }},`$ (19)
where

$`\sigma =\left({\displaystyle \frac{2}{3}}\right)^{1/2}{\displaystyle _0^x}{\displaystyle \frac{dx^{}}{\varphi (x^{})}},`$ (20)
and $`\varphi (x)`$ is the normalized Voigt function given by

$`\varphi (x)={\displaystyle \frac{1}{\sqrt{\pi }}}H(x,a).`$ (21)
Here the boundary conditions are

$`J(\tau ,\pm \mathrm{})=0`$ (22)
$`J(\pm \tau _0,x)=\pm 2H(\pm \tau ,x)={\displaystyle \frac{2}{3\varphi (x)}}\left({\displaystyle \frac{J}{\tau }}\right)_{\pm \tau _0}.`$ (23)
Since the source function $`E(\tau ,x)`$ in the right hand side of Eq. (19) is very sharply peaked, we can apply the Green function method to the problem. With an approximaton

$$\varphi (x)a/\pi x^2$$
(24)
and an application of the Sturm-Lioville theory, Neufeld (1990) derived a solution for the case of a mid-plane source,

$`J(\pm \tau _0,x)={\displaystyle \frac{\sqrt{6}}{24}}{\displaystyle \frac{x^2}{a\tau _0}}{\displaystyle \frac{1}{\mathrm{cosh}[(\pi ^4/54)^{1/2}(|x^3x_i^3|/a\tau _0)]}}.`$ (25)
We perform Monte Carlo calculations for a monochromatic source with $`x_i=0`$, and compare the result with Eq.(24). Fig. 4 shows our results for $`\tau =10^5`$ and $`\tau _0=10^6`$. The solid lines are the results of Monte Carlo calculation, and the dotted lines are those obtained by Neufeld. The shapes of the profiles agree with each other, and we emphasize that our Monte Carlo code is written to incorporate all the quantum mechanics associated with both resonant and non-resonant Ly$`\alpha `$ scattering.

The figure shows an overall disagreement: the Monte Carlo solution appears to be translated from the Neufeld’s solution by an amount of $`x3`$, which corresponds to the frequency shift where the wing scattering becomes important and the profile function $`\varphi (x)`$ can be nicely approximated by Eq. (24). In fact, in the cases of dust free media, no loss of Ly$`\alpha `$ line photons is permitted, which guarantees the flux conservation. The profile function $`\varphi (x)`$ can not be approximated by Eq.(23) at $`x3`$, and therefore, Neufeld’s calculation underestimates the amount of core photons that are removed and ultimately redistributed to the wing regimes, which reduces the amount of diffusively trasferred wing photons.

### 3.3 Survival Fraction for Dusty Media

Neufeld (1990) also provided an analytic solution for the problem of the radiative transfer in a dusty medium. Here we will concentrate on the continuum absorption case. Harrington (1973) and Neufeld (1990) formulated this problem by

$`{\displaystyle \frac{^2J}{\tau ^2}}+{\displaystyle \frac{^2J}{\sigma ^2}}=3[\varphi (x)+\beta _t(x)]\left[{\displaystyle \frac{E(\tau ,x)}{4\pi }}\beta (x)J\right],`$ (26)
where

$`\sigma =\left({\displaystyle \frac{2}{3}}\right)^{1/2}{\displaystyle _0^x}{\displaystyle \frac{dx^{}}{\varphi (x^{})+\beta _t(x^{^{}})}}.`$ (27)
We assume that the total dust opacity $`\beta _t(x)0`$ becasue $`\varphi (x)\beta _t`$. Here, $`\beta `$ is an absorptive part of dust opacity per atomic hydrogen, which is given by

$`\beta =1.0\times 10^8T_{g4}^{1/2}{\displaystyle \frac{n_H}{n_{HI}}}\xi _d.`$ (28)
Here

$`\xi _d=(1A){\displaystyle \frac{\tau _d}{N_{HI}}}/\left[(1A){\displaystyle \frac{\tau _d}{N_{HI}}}\right]_G,`$ (29)
where A is albedo and the subscript $`G`$ represents the Galactic value. A dust opacity per a hydrogen atom is given by

$`{\displaystyle \frac{\tau _d}{N_{HI}}}=f_d(\sigma _{sca}+\sigma _{abs}),`$ (30)
where $`\tau _d`$ is dust optical depth, $`f_d`$ is dust-to-gas ratio, $`\sigma _{sca}`$ and $`\sigma _{abs}`$ are the scattering and the absorption cross sections, respectively. According to Draine and Lee (1984), the mean value of the dust optical depth for Our Galaxy is given by $`\tau _d/N_{HI}=1.6\times 10^{21}`$ at $`\lambda _0=1216\AA `$.

In our Monte Carlo code, we give $`\tau _d`$, $`A`$, as well as $`N_{HI}`$ or $`\tau _0`$. Therefore we can compare these two kinds of calculations. Neufeld also gave solutions for two kinds of approximation: (1) variable effective scattering frequency, (2) constant effective scattering frequency. In the former approximation, the survival fraction of photons is given by

$`f_e=\left[1+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{y_0^{2k}}{k!}}{\displaystyle \frac{1}{3711\mathrm{}(4k1)}}\right)\right]^1,`$ (31)
where $`y_0\frac{3}{2}(\frac{16}{\pi })^{1/6}a^{1/6}\beta ^{1/2}|\tau _0|^{2/3}`$.

In the latter approximation, the total fraction of photons which escape the slab is given by

$`f_e={\displaystyle \frac{1}{\mathrm{cosh}(Y_0)}},`$ (32)
where

$`Y_0=\left[{\displaystyle \frac{3(a\tau _0)^{1/3}\beta \tau _0}{\pi \zeta ^2}}\right]^{1/2},`$ (33)
and $`\zeta `$ is determined to be $`\zeta =0.525`$ by fitting the results of Hummer and Kunasz (1980) to the following relation between peak frequency ($`x_s`$) and the $`a\tau _0`$,

$`x_s=\zeta (a\tau _0)^{1/3}.`$ (34)
The other way of determining $`Y_0`$ is provided by Neufeld (see Eq. 4.35 in his paper), who considered the mean scattering length in the result of Hummer and Kunasz (1980). The consideration leads to an survival fraction $`1/\mathrm{cosh}(Y_0)`$ with

$`Y_0=\left(31.8\beta \tau _0\right)^{1/2}.`$ (35)
He argued that this relationship may be applicable even in the range of intermediate optical depth, $`10^3<\tau _0<10^3/a`$.

In Fig. 5 we show the main result on the survival fraction of Ly$`\alpha `$ photons in a dusty medium with $`\tau _0=10^6`$ and $`T=10^4\mathrm{K}`$. Our results are in good agreement with Neufeld’s approximate solutions and those of Hummer and Kunasz. Our results are also consistent with Neufeld’s approximation solution for the approximation of a constant effective frequency.

On the other hand, although Neufeld noted that Eq.(34) above is valid even in intermediate optical depths, according to our results, this statement is not true.

## 4 Discussion

In this paper we described the Ly$`\alpha `$ line transfer in an optically thick, dusty, and static medium. We took accurate atomic physical considerations in the Monte Carlo code, and examined the effect of dusts on Ly$`\alpha `$ line formation in optically thick media.

For dust free media we confirmed the line transfer mechanism, in which a number of longest flights occur. We emphasize that photons experience a series of wing scatterings at the moment of escape, at which stage significant polarization can be developed. We also confirmed that line profiles from our Monte Carlo calculations agree well with the analytic solution derived by Neufeld (1990).

For dusty media, we calculated the survival fraction of photons for $`\tau _0=10^6`$ and $`T=10^4K`$, and found that our results are in good agreement with those obtained by Hummer and Kunasz (1980) and with those of Neufeld’s approximation solution.

The computational speed is limited for extremely thick media with $`a\tau _0>10^3`$. This difficulty can be overcome by considering only wing scatterings (Adams 1972). Hence, we are now writing a Monte Carlo code similar to that of Adams. Also we will investigate effects of following various deviations from our simplified configuration: the initial line width, the bulk motion of scatterers, the geometrical shape of scattering media, the degree of homogeneity of the dusts, and the distribution form of photon sources.

One of main goals of our research is to develop a tool to estimate the dust abundance in accordance with the usual method using the Balmer decrement (see Osterbrock 1989) and the slope of UV continuum (Fall & Pei 1989). This may help us correct properly the dust extinction of emission lines in starburst galaxies. Especially this will enable us to estimate the star formation rate of those galaxies which are located near or far.

The other ramification is the formation of the damped Ly$`\alpha `$ absorption (DLA) in many quasar spectra. The dust contents in DLA galaxies may deform the absorption profile. This may spoil the usual Voigt profile fitting procedure, and subsequent derivation of physical values can be erroneous.

HML and SHA are grateful for the financial support of Brain Korea 21 project. HML was also supported by KOSEF grant no. 1999-2-113-001-5. HWL gratefully acknowledges support from the BK21 project.


format: text

text:
# SOFT PARTON RESUMMATION IN THE CURRENT REGION OF SEMI-INCLUSIVE DEEP INELASTIC SCATTERING††footnote †Presented at the 8th Intl. Workshop on Deep Inelastic Scattering (DIS2000), Liverpool, U.K., April 2000

## Acknowledgments

P. N. thanks the organizers of the DIS2000 Workshop for their warm hospitality. This work was supported in part by the NSF under grant PHY-9802564.

## References


text:
# Measuring Supernova Neutrino Temperatures using Lead Perchlorate

## I Introduction

Recently a number of groups have expressed interest in using Pb as a target for neutrino interactions to study supernovae or oscillations, and this inspired an estimate of the cross section. The interest arises because of the large cross section and the low relative cost of Pb. As a result, additional cross section calculations were done recently by Fuller, Haxton, and McLaughlin (hereafter referred to as FHM) and Kolbe and Langanke (hereafter referred to as KL).

The interesting neutrino interactions with Pb consist of

$`\begin{array}{cccc}\nu _e+^{208}Pb\hfill & & {}_{}{}^{208}Bi^{}+e^{}\hfill & \hfill (CC)\\ & & \hfill & \\ & & {}_{}{}^{208y}Bi+x\gamma +yn\hfill & \\ & & & \\ \nu _x+^{208}Pb\hfill & & {}_{}{}^{208}Pb^{}+\nu _x^{^{}}\hfill & \hfill (NC)\\ & & \hfill & \\ & & {}_{}{}^{208y}Pb+x\gamma +yn.\hfill & \end{array}`$

Fig. 1 shows the energetics of these transitions. The number of neutrons emitted (0, 1, or 2) depends on the neutrino energy and whether the transition is induced by a charged current (CC) or neutral current (NC) interaction. The nuclear physics of this system is described in detail in Refs. .

An ideal lead-based neutrino detector would have an appreciable density of lead atoms and the capability of detecting the electrons, gammas and neutrons produced in the reaction. Lead perchlorate (Pb(ClO4)2) has a very high solubility in water (500 g Pb(ClO4)2/ 100 g H2O ) and the mixture is transparent. This transparency raises the hope that a $`\stackrel{ˇ}{C}`$erenkov detector can be assembled. Additionally, the presence of a neutron moderator (hydrogen) and a neutron capture nucleus (35Cl) provides a method for observing the neutrons. The high-energy (8.6 MeV), neutron-capture $`\gamma `$ rays from Cl would Compton scatter in the fluid and be observed via the $`\stackrel{ˇ}{C}`$erenkov light of the recoil electrons.

Supernovas emit on the order 1058 $`\nu `$ of all types. The average energies of neutrinos emitted by a supernova follow a hierarchy: $`\overline{E}_{\nu _e}<\overline{E}_{\overline{\nu }_e}<\overline{E}_{\nu _{\mu ,\tau }}`$ where $`\nu _{\mu ,\tau }`$ indicates either $`\mu `$ or $`\tau `$ neutrinos and their anti-particles. The predicted values of the average energies for the three neutrino classes differ a great deal and are estimated to be $`\overline{E}_{\nu _e}=11`$ MeV, $`\overline{E}_{\overline{\nu }_e}=16`$ MeV, and $`\overline{E}_{\nu _{\mu ,\tau }}=25`$ MeV.

Although not strictly thermal, the supernova neutrino energy spectra can be described well by a thermal distribution,

$`f_\nu `$ $`=`$ $`{\displaystyle \frac{1}{T_\nu ^3F_2(\eta )}}{\displaystyle \frac{E_\nu ^2}{exp(E_\nu /T_\nu \eta )+1}}`$
where $`f_\nu `$ is the normalized neutrino spectral shape, $`T_\nu `$ is an effective neutrino temperature, $`E_\nu `$ is the neutrino energy, and $`\eta `$ is the degeneracy parameter (chemical potential divided by $`T_\nu `$). For $`\nu _e`$ and $`\overline{\nu }_e`$, the supernova model predictions for the flux is well fit with $`\eta `$ = 3, whereas for $`\nu _{\mu ,\tau }`$ the flux is better described with $`\eta `$ = 0. The normalization factor $`F_2(\eta )`$ has the values $`F_2(0)`$ = 1.803 and $`F_2(3)`$ = 18.969. The total neutrino fluence $`F_\nu `$ from a supernova can be estimated by

$`F_\nu `$ $`=`$ $`2.8\times 10^{11}cm^2{\displaystyle \frac{E_B}{10^{53}erg}}{\displaystyle \frac{MeV}{T_\nu }}({\displaystyle \frac{10kpc}{D}})^2`$
where $`E_B`$ is the released energy and D is the distance to the supernova.

The supernova $`\nu _e`$ spectrum is too soft to produce a large number of multiple-neutron inverse-beta-decay events in Pb. However, if the higher energy $`\nu _{\mu ,\tau }`$ oscillate into $`\nu _e`$, the hardening of the spectrum will greatly alter the expected number of 2-neutron events. Thus the response of a Pb detector is intriguing because if T$`_{\nu _e}`$ is found to be larger than T$`_{\overline{\nu }_e}`$, it would be strong evidence that neutrino oscillations are taking place. However, the uncertainties in the neutrino energy distribution as predicted by the supernova models may complicate the interpretation of a measurement of the spectral parameters.

The FHM paper demonstrated that the NC and CC interaction rates in Pb are sensitive probes of $`T_{\nu _\mu ,\nu _\tau }`$, $`T_{\nu _e}`$, and $`T_{\overline{\nu }_e}`$. This paper further discusses the details of studying supernova neutrinos with Pb. In particular the use of the added data arising from observing the product electrons and gammas in coincidence with neutrons is described. These data permit the separation of the NC and CC rates and a further division of CC $`\nu _e`$ and $`\overline{\nu }_e`$ events. This collection of data can be interpreted in terms of the 3 temperatures of interest. The CC event division is discussed here in detail. Although any Pb-based experiment that measures electrons, gammas and neutrons might apply the techniques described herein, the example of using lead perchlorate as a target material is discussed.

## II Physical and Optical Properties

To build a reasonably large detector viewed by photo-multiplier tubes from the periphery, the attenuation of the $`\stackrel{ˇ}{C}`$erenkov light must be minimal. Therefore the transmission of light through the medium is critical. To understand the transmission, measurements were made at several wavelengths through several concentrations of Pb(ClO4)2. The maximum achievable transmission is still under study with attentuation lengths exceeding 2.5 m having been measured in an 80% solution. A large mass detector may require segmentation to minimize adverse effects from light loss.

The $`\stackrel{ˇ}{C}`$erenkov response of the medium depends on the index of refraction, which for Pb(ClO4)2 is 1.5 for an 80% solution by weight. The density of this solution is 2.7 g/cm3, the corresponding 208Pb number density is $`1.7\times 10^{21}`$/cm3 and the hydrogen number density is $`3.6\times 10^{22}`$/cm3. Some additional pertinent data for Pb(ClO4)2 are given in Table I.

The number of $`\stackrel{ˇ}{C}`$erenkov photons emitted by an energetic electron in Pb(ClO4)2 is about 185/cm in the wavelength region of interest for phototubes. The stopping power of 80% (50%) Pb(ClO4)2 is about 0.2 cm/MeV (0.33 cm/MeV). Hence a 15 MeV electron will emit a total of about 550 (920) photons in 80% (50%) Pb(ClO4)2. Due to the modest light yield, the resolution ($``$20%) will be modest. The resonance widths will also contribute an uncertainty in relating the electron energy to the incoming neutrino energy by a comparable amount (few MeV).

Because solutions of Pb(ClO4)2 contain a great deal of water and Cl, the neutron thermalization and capture time is of the order of 10-100 $`\mu `$s. This time scale is short compared to the $`\nu `$ detection rate during a supernova, but long enough to identify whether there are 0, 1, or 2 neutrons in coincidence with a primary e- or $`\gamma `$.

## III Cross Section Models

The FHM paper pointed out the importance of the 208Pb \- $`\nu `$ cross section and the production of neutrons. They noted that the total number of neutrons produced by neutrino interactions with Pb is sensitive to the effective temperature of the supernova $`\nu _e`$ energy spectrum especially for reactions which produce multiple neutrons. However, they did not consider detection schemes that could measure the energy released in the form of electrons and gamma rays (hereafter referred to as the electromagnetic energy). Neutral current events release a large number of neutrons but very little electromagnetic energy. Hence NC events can be isolated from CC events by measuring the electromagnetic energy in coincidence with the neutrons.

The papers by FHM and KL provided effective cross sections for neutrino-lead interactions averaged over the neutrino energy spectra and summed over all the product nucleus states. It is just these distributions, however, which are needed to simulate the response of a detector and to optimize its design parameters. This section summarizes the nuclear physics model implemented to calculate these distributions. In the following section, these distributions and their use in analyzing the neutrino spectra are described.

Specifically, the energy dependence of the five primary transitions of interest as identified by FHM is calculated: two for the NC interaction and three for the CC interaction. The parameters describing these transitions are given in Table II. For the isobaric analog state (IAS) CC and Gamow-Teller (GT) NC and CC transitions, the analytic expressions and matrix element values from Ref.  are used. For the first forbidden transitions, the cross section can not be written in closed form. Software to calculate that contribution to the cross section as used by FHM was obtained and hence the cross section calculations presented here agree with FHM by design.

Figure 2 shows the cross section overlaid on a normalized thermal flux distribution ($`T_{\nu _e}`$=2.76 MeV, $`\eta `$=3). Figure 3 shows the same for a higher neutrino temperature ($`T_{\nu _e}`$=6.27 MeV, $`\eta `$=3). The origin of the large increase in the interaction rate with temperature is clear from a comparison of these plots. The onset of the Pb cross section lies between these two temperature extremes. As a result, the interaction rate in Pb is highly sensitive to $`T_{\nu _e}`$. This is especially true for the first forbidden transition which leads to the emission of 2 neutrons.

Tables III and IV show some cross sections and relative interaction rates for the various components of Pb(ClO4)2. It is clear that the interactions with Pb and H dominate. Contributions from Cl, O, and elastic scattering of electrons are at the few percent level. The CC interaction of $`\overline{\nu }_e`$ on protons is a major contribution because the hydrogen density is large in Pb(ClO4)2 solutions. For other isotopes found in Pb(ClO4)2 but not listed in the table, FHM note that the cross section models vary like N-Z or NZ/A. For Pb, these dependencies are small for the naturally occurring isotopes. Hence one can anticipate that their response to neutrinos will be similar. This is not the case for Cl. However, the difference in cross section between 35Cl and 37Cl is still probably smaller than other uncertainties and thus the Cl interaction rate is small.

## IV The Electromagnetic Energy

A detector based on lead perchlorate would be able to measure the energy deposited by e- and $`\gamma `$ in coincidence with the number of neutrons produced. The energy deposited in the detector due to e- and $`\gamma `$ is here referred to as the electromagnetic energy to contrast it with that resulting from neutron capture. To identify a particular class of neutrinos, one sorts the electromagnetic energy spectra by the number of coincident neutrons. CC events always produce an electron which can have substantial energy in coincidence with neutrons. NC events produce either a high-energy $`\gamma `$ with no neutrons or neutrons with little or no $`\gamma `$ energy.

First consider the separation of NC from CC events. The product nucleus from the neutrino interaction can produce 0, 1, or 2 neutrons. (Particle decays other than neutrons (alphas and protons) of the product nucleus states have a much lower branching ratio (see KL) and are ignored here.) When the interaction does produce a neutron, the neutrons carry away the available energy and there is little electromagnetic energy produced. Hence for the NC events, the electromagnetic energy is only appreciable when 0 neutrons are emitted. Thus events which do have appreciable electromagnetic energy in coincidence with neutrons are indicative of 208Pb - $`\nu _e`$ or H - $`\overline{\nu }_e`$ CC interactions.

Next consider the separation of the CC events into $`\nu _e`$ or $`\overline{\nu }_e`$ events. The 2-neutron event rate is dominated by $`\nu _e`$-induced CC transitions to the first forbidden level in Bi. Since the electromagnetic energy for the CC reactions is related to the incoming neutrino energy, the 2-neutron event spectrum provides data on the incident $`\nu _e`$ energy spectrum. In contrast, the IAS and GT transitions induced by $`\nu _e`$ produce mostly 1-neutron events. Furthermore, although $`\overline{\nu }_e`$ do not have an appreciable CC cross section on Pb, they have a large cross section on H. The 1-neutron spectrum is therefore generated by reactions induced by $`\nu _e`$ and $`\overline{\nu }_e`$. But those are the only reactions that contribute significantly to the 1-neutron electromagnetic energy spectrum and therefore the 1-neutron spectrum provides data on the $`\overline{\nu }_e`$ energy distribution. The 2-neutron and 1-neutron spectra thus relate directly to the spectral features of the incident neutrino flux. Figs. 4 through 7 show spectra of the electromagnetic energy for events in coincidence with 1 or 2 neutrons for two different neutrino energy distributions. These spectra have been convolved with an estimate for the detector resolution.

For NC events, only GT excitations can produce an event with no neutrons as the resonance straddles the neutron separation energy. A $`\gamma `$ ray of approximately 7.3 MeV is then emitted. This would produce a line in the 0-neutron electromagnetic energy spectrum. NC transitions can also produce 1-neutron events and a modest number of 2-neutron events. As there will be little electromagnetic energy in coincidence, the 1-neutron events will be interpreted as an 8.6-MeV line in the 0-neutron spectrum. The strengths of these two lines in the 0-neutron spectrum provide a measure of the NC rate. Due to the resolution, these two lines will most likely be blended. These rates will be large even for low-temperature profiles.

One background to the CC 1-neutron spectrum comes from 2-neutron NC events. One of the 2 neutrons would be taken as a neutron but the second would be mistaken for electromagnetic energy. For $`T_{\nu _\mu ,\nu _\tau }`$=6.27, $`T_{\nu _e}`$ =2.76, $`\eta `$=3, the NC rate of 2-neutrons events will be comparable to the CC rate because, although the NC cross section is much smaller, there are 4 times as many $`\nu _{\mu ,\tau }`$ as $`\nu _e`$ and their temperature is much higher. Thus at low $`T_{\nu _e}`$, this background must be separated. Note that because there are several $`\gamma `$ sharing the energy and photons Compton scatter more than once before losing all their energy, the $`\stackrel{ˇ}{C}`$erenkov light produced by neutrons capturing on Cl will be isotropic in comparison to electrons from inverse beta decay. This could lead to a technique for separating this background. Otherwise an electromagnetic energy threshold of 10 MeV would eliminate it.

With the spectra in hand, one deduces the $`\nu `$ temperature by relating it to the average electromagnetic energy from the 1-neutron spectrum ($`\overline{E_1}`$) and the 2-neutron spectrum ($`\overline{E_2}`$) (Fig. 8). One also counts the total number of 1-neutron events ($`N_1`$) and 2-neutron events ($`N_2`$) in these spectra (Fig. 9). The ratio (R=$`N_2`$/$`N_1`$) depends strongly on the temperatures of the two types of neutrinos but is insensitive to uncertainties in the supernova distance scale. Using these parameters ($`\overline{E_1}`$, $`\overline{E_2}`$, $`N_2`$, $`N_1`$, R), the pair of neutrino temperatures that fit the data are determined.

If the electron neutrino temperature is clearly greater than that of the electron anti-neutrino temperature, it would be considered strong evidence for $`\nu _e\nu _\mu `$ or $`\nu _e\nu _\tau `$ oscillations. This condition is satisfied by the top left hand section of these two-temperature plots. It should be noted that the contour plots were done for thermal spectra. If $`\nu `$ oscillate, the spectral shape could be complicated depending on the oscillation conditions and this analysis would have to be extended.

To use this set of figures to determine the temperatures with some precision, one must estimate the uncertainty in the measured parameters. These uncertainties are mostly determined by the number of events observed. Regions in Figs. 8 through 10 are defined by the values of these parameters and their uncertainties. The boundaries of these regions indicated on the axes provide the corresponding temperature ranges that describe the data. For a commonly considered supernova at 10 kpc with $`3\times 10^{53}`$ ergs of released energy and a 10 kt-Pb detector (24 kt of lead perchlorate), one expects to measure $`\overline{E_2}`$ = 11.0 $`\pm `$ 0.5 MeV, $`\overline{E_1}`$ = 19.6 $`\pm `$ 0.2 MeV, $`N_1`$ = 1530, $`N_2`$ = 82, and R = 0.053 $`\pm `$ 0.006 for non oscillating spectra ($`T_{\nu _e}`$ = 2.76, $`\eta `$ = 3, $`\overline{E}_{\nu _e}`$ = 11; $`T_{\overline{\nu }_e}`$ = 4.01, $`\eta `$ = 3, $`\overline{E}_{\overline{\nu }_e}`$ = 16). The uncertainties for the average energies were estimated by calculating the variance of the distribution and dividing by the square root of the number of observed counts. The resulting statistical uncertainty in the temperatures is small and the precision improves if $`T_{\nu _e}`$ is larger because the event rate is higher.

## V Summary

The Sudbury Neutrino Observatory (SNO) and Super-Kamiokande experiments are operating. In the event of a supernova, Super-Kamiokande will observe a large number of $`\overline{\nu }_e`$-H events ($`10000`$ for $`T_{\overline{\nu }_e}`$ = 4.01, $`\eta `$ = 3) and SNO will see a number of $`\nu _e`$-D events ($`90`$ for $`T_{\nu _e}`$ = 2.76, $`\eta `$ = 3). The average energy measured for these two reactions could be compared in a two-temperature plot in a similar analysis as described above. Due to the smaller number of events, the uncertainty comparing these energy distributions would be determined by the SNO measurement. For these $`\nu `$ temperatures, a 24 kt Pb(ClO4)2 experiment would count a number of 2-neutron $`\nu _e`$ events comparable to the number of SNO D($`\nu _e`$,e-)pp events. Thus the average energy analysis of a Pb(ClO4)2 experiment would have a sensitivity similar to the SNO-SK comparison. But the total number of counts in a Pb(ClO4)2 experiment is much more sensitive to the neutrino temperatures. SNO (Pb(ClO4)2 ) would observe about 3 (60) times as many D($`\nu _e`$,e-)pp (2-neutron) events at $`T_{\nu _e}`$ = 6 than at 2.76 MeV. This extra handle on the data adds a powerful redundancy to the analysis.

It is obvious from Table V that there are still large uncertainties in the nuclear physics of these interactions. Furthermore, reactions on the other Pb isotopes have not been considered although they would tend to increase the signal. Therefore the response of such a detector should be measured in a prototype. The $`\nu _e`$ spectrum from stopped $`\mu ^+`$ decay is comparable to that from supernovae. So the 2-neutron and 1-neutron spectra could be measured at a beam stop facility  without the influence of $`\overline{\nu }_e`$. Since the $`\overline{\nu }_e`$-proton reaction is well understood, it does not need additional study and the measurement would focus on what is currently poorly known: the nuclear physics parameters of the lead reactions. Thus the theoretical uncertainties in the cross section could be diminished by a beam stop experiment.

The details of the analysis presented in this paper will alter as the nuclear physics of this system becomes better known. But the general conclusions generic to lead detectors will not. It is very powerful to measure not just the total number of neutrons produced but the 2-neutron and 1-neutron rates separately. Furthermore the observation of electrons and gamma rays in coincidence with the neutrons can permit the separation of NC and CC interactions. If one observes a supernova with a Pb detector, these combinations of data can be analyzed to estimate the temperatures of the neutrino spectra.

## VI Acknowledgements

I would like to thank Peter Doe, Hamish Robertson, Thomas Steiger, John Beacom, Edwin Kolbe, Gail McLaughlin and George Fuller for useful conversations. I am especially grateful to Wick Haxton for providing code to estimate the first forbidden transition cross section. This research was support by the Department of Energy, grant number DE-FG03-97ER41020.


text:
# A Novel Method for Sampling Alpha-Helical Protein Backbones

## I Introduction

### A Outline of current prediction methodology

Prediction of protein structure from sequence is one of the most enticing goals of scientific inquiry today. Currently the most reliable method of determining a protein’s shape is to search for a close homologue in the database of solved protein structures. Although the number of solved structures increases daily, it is estimated (Brenner et al., 1997) that in the near future at least 40% of proteins of interest bear no discernible sequence resemblance to a known macromolecule. Therefore ab initio prediction of structure from sequence remains an important challenge.

When sequence homology cannot be used to construct a 3-D model, the current modus operandi for predicting structure is composed of three separate (yet interdependent) steps: sampling, searching, and ranking. First one picks a formalism to index and sample the possible structures, then searches through the conformation space and constructs the shapes, and finally uses some ranking criterion to pick out a structure as close to native as possible. Let us provide a few examples of each step. Some of the possible choices in the initial phase, picking a representation, are all-atom, reduced ’pseudo’ amino-acid interaction centers (Levitt & Warshel, 1975; Levitt, 1976), lattice (Hinds & Levitt, 1992; Covell, 1994; Lau & Dill, 1989; Skolnick & Kolinski, 1990), and many others. Once the representation is chosen, one then picks the degrees of freedom to use. Researchers may vary Cartesian coordinates of all or a subset of the atoms, dihedral angles, relative distances, etc (Eyurich et al., 1999). The decisions of which, how many, discretized or continuous, and so on, are motivated by the amount of detail one wants to capture, as well as by computational complexity of the forthcoming search.

Having picked a representation, one then selects an appropriate search technique, or a mixture of search techniques to march through the enormous space of possible 3-dimensional structures. As we have mentioned above, the three complementary steps of prediction are highly interdependent. The most important factors in picking a search method are computational complexity and the natures of both the energy function and the representation. Lattice models, for instance, lend themselves easily to exhaustive enumeration. All atom Cartesian coordinate representations usually employ a physical potential, and search for a minimum with either an adaptive minimization routine, or a copy of Nature’s algorithm, molecular dynamics.

The remaining tool employed in protein structure prediction is a potential function constructed to pick out native or near-native conformations from a vast number of alternatives. If the search method employed is a guided one (for example, minimization) then the potential function also provides a landscape that will allow the procedure to converge to the correct answer. Just like sampling and search techniques, potential functions are many and varied: physical potentials with/without water representation, knowledge based potentials of various detail, hydrophobic contact potentials, and many others (Sippl, 1990; Samudrala, 1997; Park & Levitt, 1996). Frequently more than one potential is used to evaluate candidate structures.

## II Description of the Sampling Method

### A Graph-Theoretical representation of secondary structure

We represent alpha-helical proteins as connected graphs. Each helix is represented by a vertex and an edge is drawn between two vertices if the corresponding helices are in contact. To form a protein core one builds a connected subgraph, adding one helix at a time, until the protein is assembled. Figure 1 shows the graph-theoretical representation and the construction of Myoglobin.

We are, in essence, building an off-lattice model for secondary structure segments, but our degrees of freedom are not in the frequently chosen (Park & Levitt, 1996; Eyrich et al., 1999) $`(\varphi ,\psi )`$ loop residues (or their subsets), but in the relative geometric position of the helices themselves. The disadvantage of this approach is that a significant fraction of visited structures violate loop constraints - simply put, the ends of helices are too far apart to be joined together by an intervening loop. There are, however, several advantages. First, the structures we generate have a much higher tendency to be compact. Second, we are able to exploit the correlation between the sequence patterns of helix-helix contact to significantly enrich our sampling with native-like structures. (A simple example of this is two helices with small residues in the contact area are more likely to be close to each other than those with large residues between them.) Third, we are able to sample possible packing much better than when using loop torsion angles, as is usually done in off-lattice models.

### B Exhaustive enumeration of helix-helix contacts

In this work our principal aim is to achieve coarse-grained sampling of the helical protein core. Our hypothesis is that the contact patches on each helix influence the way the helices pack. We choose to derive the preferences for helix-helix packing from the distribution of such orientations in known proteins. If we are incorrect, then our packing strategy will position helices randomly and no harm will be done. Therefore, the way we choose to travel through our space is by exhaustive enumeration of a discrete version of our representation. We can choose to sample relative helix-helix orientation in several ways. (We call a particular realization of an edge a ’link’. Figure 2 shows the various attributes of a link). It’s useful, once again, to draw an analogy to dihedral sampling: there one can discretize a particular dihedral angle into a (not necessarily) uniform spectrum of values suggested by geometric considerations; e.g. sampling a dihedral angle in, say, 30 degree intervals. Alternatively, one can extract the local moves from a database of existing structures, thus sampling the space more efficiently. One could, for example, assign three possible values - helix, loop, sheet - to a dihedral angle. Another possibility would be to bias the assignment with sequence matching. (Simons et al. 1997) or secondary structure prediction. In this work we choose exhaustive enumeration of possible orientations derived from known structures. This is motivated by the observed influence sequence has on the packing of proteins (Reddy & Blundell, 1993; Popov, 1980; Matheson & Sheraga, 1978; Richmond & Richards, 1978).

### C Definition of helix to helix contact (’link’) and the contacts’ database

We extract the possible relative orientations of two helices from a database of touching helix-helix pairs obtained from a subset of the SCOP (Murzin et al.,1995) database. The 1305 folds in our library have a sequence identity to each other no greater than 35%. We then parse our database to get touching helix-helix pairs. We define two helices $`X`$ and $`Y`$ to be in contact if a) the shortest distance between any two CB’s, $`CB_x`$ and $`CB_y`$, located respectively on helix $`X`$ and helix $`Y`$, is less than 7.3 Å; and b) if each helix has at least 3 CB atoms within $`d=d_{min}+2.5`$ Å. ***The three CB on each helix are necessary to define a relative rotation from one helix to another. The specific values of $`d_{min}=7.3`$Å and 2.5Å worked best in our tests. Lowering the $`\mathrm{d}_{min}`$ parameter to below 5Å caused some of the proteins in our database to be represented by disconnected graphs, thus making it impossible to ever reproduce them with our technique. Conversely, making the $`\mathrm{d}_min`$ cutoff much larger than 7.3$`\AA `$ decreases the influence the contact sequence has on the relative orientation of the helices. In addition to the minimal three CB’s on each helix, all the CB’s on both helices that are within 2.5 Å of $`d_{min}=7.3`$ are defined to be in a contact ’patch’. Our definition picks out the residues in the contact region by assigning a ’patch’ of contact residues between helices. Residues on the far side of each helix, which have very little influence on the relative orientation, are ignored. Figure 2 shows the patches and the structure of each link in the database.

### D The Enumeration Procedure

Figure 1 illustrates the enumeration procedure. The build-up of a helical protein proceeds as follows: given a protein sequence and native secondary structure assignment, we construct idealized helices on the chain. (In a predictive scenario the native secondary structure assignment will be replaced by a prediction, or possibly several alternative predictions.) We then go to the library of links and perform a sequence alignment between the patches on each link and all possible pairs of helices on our target sequence. The residues between key ’patch’ residues serve as spacers to fix the position of the influential patch residues. If the sequence match score is high then this particular link will be used to bring the pair of helices together. The actual threshold for matching depends on the number of structures we want to sample. The sequence match is scored using the Blosum62 (Henikoff & Henikoff, 1992) matrix. The comparison threshold typically varies from 0.1 to 0.6 and is adjusted to give anywhere from 10 to 1000 possible links for a given pair of helices in the target protein. The number chosen depends on how many final structures we wish to generate in the available amount of computer time. The relative orientations of the helices are then loaded into memory and the buildup of structures begins.
We use each topologically distinct pathway to build the target - (refer to figure 1 b) ) - each specific pathway corresponds to a minimal subgraph spanning the protein graph. (Dashed and solid lines, respectively, in figure 1). For a structure consisting of 4 helices this results in $`4^{(42)}=16`$ possible topologies (see Figure 3). Note that the graph corresponding to each final structure might (indeed should) possess other edges not used in construction; however it’s sufficient to follow a minimally connected subgraph to construct it.

For each topology we construct each possible combination of links (gathered from the sequence-matching procedure above) that can realize a given edge. To improve performance we used branch and cut filters for loop and clash constraints. The clash filter eliminates a conformation if more than 3 residues on one helix are closer than 2 $`\AA `$ to residues on another helix. The loop filter eliminates conformations for which the distance for the loops necessary to connect the helices is longer than the maximum available loop length. The reason for the branch-cut approach is simple: if in a given six-helix enumeration helix one clashes with helix three there is no need to cycle through and build helices four, five, and six. Finally, each geometrically viable structure is tested for compactness. All the tests and filters are extremely fast because whenever possible we use the coarse segment representation of the structure and thus escape having to visit each amino acid’s coordinates. We have also incorporated other filtering information, most notably disulfide bond locations, into the build-up procedure. In its current incarnation the method is able to generate roughly $`10^3`$ conformations per second on a 400MHz Pentium workstation. The ultimate speed of the procedure will, in the future, be limited by scoring function evaluations.

## III Results

To test the performance of our technique we have used the coordinates of helical cores for 25 proteins. The molecular sizes range from 31 to 172 amino acids, and the number of helical residues to which coordinates are assigned ranged from 23 to 130.

The position of disulfide bonds can sometimes be easily obtained through chemical methods, and one can then rely on knowing their location prior to prediction of structure. To this end, for 4 of the 25 proteins in our set that have disulfide links, we have also evaluated the performance of our method with and without a priori knowledge of these bonds.

The fact that we assign residues only to the helical backbone leaves open the possibility that our results will be degraded by subsequent assignment of coordinates to all remaining atoms. We investigated this effect by building all atomic coordinates for 5 of our proteins using the program SegMod (Levitt, 1992). Prior to reconstructing full coordinates we pruned the helical cores with a slightly modified Sippl-like function (Sippl, 1990), leaving 500 to 4500 structures for each protein. (We left more decoys for larger proteins).

### A Sampling of helical cores

We used the method to generate conformations for 25 helical proteins, ranging in size from 31 to 172 residues and containing anywhere from two to six helices. The program STRIDE (Frishman & Argos, 1995) was used to determine the secondary structure of the native protein. To check how much the quality of our sampling will be degraded when all the coordinates are reconstructed from the helical core, we have included several proteins with significant loop content. The assignments of helices to the structure were made identical to the native structure.

The results are summarized in table I. Overall the results are very promising. For nearly every protein a sizeable proportion of the sampled structures is within 3Å of the native. For 9 proteins the best structure produced is closer than 1Å CA RMSD to the actual structure, and is virtually identical to the protein itself. Our procedure clearly samples well enough for a suitable potential to make a successful prediction.

Our method is not only very fast, it is also very efficient, in the sense of being able to generate a large percentage of native-like conformations. We can access the efficiency by estimating how many random protein-like compact structures are needed to obtain the RMSD of our best structure, and compare that number to the number of conformations we have visited. We use the estimate from (Reva et al., 1998)

$$\frac{1}{\left(\sigma \sqrt{2\pi }\right)}_{\mathrm{}}^R\mathrm{exp}\left[\frac{\left(x<R>\right)^2}{2\sigma ^2}\right]𝑑x$$
(1)
Where, following Reva et. al., 1998, we set $`\sigma =2.0`$ and $`<R>=3.333N^{1/3}`$, where $`N`$ is the number of residues in the protein core that we have assigned coordinates to. The proximity of our best structures to native exposes the weakness in estimating probability of low RMSDs, namely that the Gaussian distribution cannot be used to effectively describe sets of conformations which are very close to native. Having said that, we still feel that equation 1 gives a good estimate of efficiency. The last column in table I shows the ratio of the number of structures we sample to the number of random structures given by formula 1.

The results from this column in table I show that our method is efficient - the sampling is strongly biased towards the native conformation. The worst performance is shown in protein 1res, where we would have done slightly better if we chose relative helical conformations at random. However for most cases we are besting the random sampling by factors of $`10^4`$ to $`10^6`$, with efficiency actually improving for large proteins - which is where it matters most.

### B Features of the sampled conformations

Figure 4 shows a typical distribution of the sampled conformations. One particularly interesting feature is the non-Gaussian tail extending towards the native structure. This is most likely the result of using sequence matching to select the possible orientations of the helix-helix pairs. The use of sequence information enriches the set with native-like structures and makes the curve decidedly non-Gaussian.

### C Using disulfide bond information

To test how well the sampling is aided by supplementary information, we took proteins with disulfide bonds and subjected the sampling for those structures to an additional filter. We did not screen 1cc5 because although it does have disulfide bonds, they occur in the loop regions and do not help filter the helical cores. Our filters rejected conformations in which the CYS CB-CB distance was greater than 8Å. This is not a very strict filter; it simply weeds out grossly incorrect topologies. The results both with and without (The results without disulfide bond information are the same as in table I) disulfide bond information are displayed in table II. The number of structures was cut by approximately a factor of 8. While the filtering did not significantly improve the average RMSD of the ensemble, the proportion of ’good’ structures was raised approximately fivefold.

### D Completing the structure

The structures described in the previous two sections possess only backbone coordinates for residues in the helical core. It is quite possible that the successes of the sampling can be washed out when coordinates are assigned to the remaining residues and side-chain atoms. To check how the quality of the sampling changes during completion we ran our sets of conformations through the program SegMod (Levitt, 1992). We did not construct full coordinates for all twenty five proteins because we simply wanted to show that reconstruction does not significantly degrade the quality of the sampling procedure. To make the test credible, we selected proteins of varying size, and with long unassigned loops. The results are shown in table III.
Once again, the resulting ensembles of decoys contain a sizeable proportion of near-native decoys. Even though we did not produce any decoys closer than 6$`\AA `$ for 2fha, we feel that the sampling of such a large (according to ab initio prediction standards) protein - 172 residues - was successful. The RMSD of the best decoy, 7.3$`\AA `$ has log-odds of 23234 (computed from last column of table III. The total ensemble for 2fha has 0.25results are good for decoy ensembles. Whether they are good enough for prediction remains to be seen.

## IV Discussion

### A Desirable features of search methods

In this paper we presented a novel sampling technique for helical proteins. To provide a context in which one can judge the merits of our procedure, we shall review the features that make a sampling method effective.

Our assertion is that for a search method to be useful in ab-initio structure prediction, the proportion of sampled conformations that are native-like should be statistically significant. Admittedly, the definitions of both “statistically significant” and “native-like” are imprecise - but clearly, the more and the closer, the better.

The need to get close to the native fold is motivated by several considerations. First, we wish to furnish good, i.e. close to native, predictions, Since we cannot select that which we do not have, having at least one native-like conformation is an obvious necessity. Second, an energy functions can only be effective in a slice of conformational space close to the actual protein . Structures can be ’wrong’ in a vast number of ways, many of which are low scoring alternatives to the native fold in a very large space. Only a small subspace close to the actual structure can be expected to exhibit the properties of an energy funnel; i.e. a region where the energy decreases uniformly as one approaches the native, and where a majority of conformations are more favorable than all the others in the available space. We can demand that an energy function correctly discriminates a 1 Å fold from all 10 Å folds, but if our best candidate model is 8 Å cRMS from native, we cannot expect an energy function to pick this over one that is 10 Å away. The need to have more than one ’good’ structure; indeed to enrich the library as much as possible with near-native choices, is brought about by the imperfections of the available selection methods. A perfect energy function should be able to discriminate a native structure from all other alternatives - this, in fact, is what nature is able to do with striking consistency.

The current generation of potentials (Samudrala & Moult, 1997; Sippl, 1990, Simons et al., 1999) often misidentifies incorrect folds, giving them energies lower than that of the native fold. By increasing the number of ’correct’ answers we increase the statistical likelihood of making a successful prediction. In addition, some current techniques use consensus information - selecting a subset and building a consensus model from it (Huang et al., 1999). These methods clearly need a sizeable native-like population to work. To give today’s energy functions all the help they can certainly use, it’s advantageous to have as many native-like folds as possible in one’s decoy set.

To make the above discussion more concrete, we should quantify the definition of being ’near-native. How close is close enough? One possible answer comes from comparing one’s technique to random sampling. Reva et. al. (1998) suggest a ’native-like’ target value of about 6 Å cRMS deviation relative to the actual structure. Other authors (Park & Levitt, 1996; Simons et al., 1997). also propose a similar distance of 4-6 Å. This number depends on the size of the protein.

### B Comparison with other search methods

Our sampling method has both merits and drawbacks when compared with other search techniques. Many classical approaches of assembling alpha-helices into a core have been proposed (Soloviev & Kolchanov, 1981; Cohen et al., 1979, Cohen & Sternberg, 1980). One alternate approach to sampling is to assemble helices using distance geometry by either embedding distance space in three dimensions (Aszodi et al., 1995; Havel, 1991), or by minimizing against restraints (Cohen & Sternberg, 1980; Kuntzet al., 1982; Richmond & Richards, 1978; Huang et al., 1999; Chelvanayagam et al., 1998; Lund et al., 1996; Mumenthaler & Braun, 1995; Ortiz et al., 1998; Smith-Brown et al., 1998). In a recent work (Huang et al., 1999) this method produced decoy libraries with the log-odds of producing the lowest structure in the -4 to -6 range Log-odds are $`1/N_{best}`$, where $`N_{best}`$ is the number of random structures needed to find at least one structure with RMSD equal to our best structure. $`N_{best}`$ is the first number in the last column of table I.. Our results, which range from log-odds of -4 to -9, are somewhat better. However we must emphasize that the comparison is unfair, since we used actual structure, and Huang et. al. used predicted (albeit well predicted) secondary structure. We are currently examining how well our method will perform with predicted secondary structure. Preliminary indications are that both distance geometry and our method both get sufficiently close to the native fold, with our method having a slight edge due to its speed.

Our method begins to pull away from distance geometry methods in the scaling of CPU time needed when one moves to larger and larger proteins. Current distance geometry methods cannot be extended to proteins of length over 100 residues (Huang et al., 1999). In contrast, because our sampling technique scales with the number of secondary structure segments and not the number of residues, it can easily handle chains of 100 - 200 residues, comparable to single domains of larger proteins. The library for the largest protein in our work - 2fha, 172 residues long - took approximately 15 hours of CPU time on a 400 MHz Pentium II machine. The sampling procedure produces approximately 500 structures per second for two-helical proteins, 100 per second for proteins containing three helices, 10 per second with four, and roughly one structure/second for proteins with five and six helices.

Another sampling technique is to fix the helical segments and then to vary the dihedral angles of the loop regions. This method has been used to produce fold libraries (Park & Levitt, 1996), and, combined with a branch and bound algorithm, in folding studies (Eyrich, et al., 1999). This method carries an advantage over ours because every structure generated automatically satisfies loop constraints, where as many of our potential structures do not. Many of our structures have to be rejected because the physical distance between ends of helices violates chain connectivity. Furthermore, we do not have loops at the end of construction; therefore we are currently working on an algorithm to sample loop regions of our protein cores that would be a good (i.e. fast) match to our helical core sampling.

Our method out-paces the loop dihedral angle sampling in two categories: scaling of computational demands with time, and efficient sampling of structures with plausible contacts. Let us illustrate the first difference with a four-helical protein. Assuming an average loop length of seven residues varying loop conformation demands 7(residues)\*2($`\varphi ,\psi `$)\*3 = 42 degrees of freedom. Our procedure needs to position three rigid bodies (the position of the initial helix is arbitrary), requiring 6(rigid body)\*3 = 18 degrees of freedom. The second difference stems from the fact each helix in our construction method is guaranteed to have at least one, and possibly two plausible contacts with other helices. When one varies dihedral angles, most of the non-clashing structures have large voids between helices.

Yet another method for generating possible reduced-model conformations for protein structure are lattice and off-lattice models, such as Hinds & Levitt, 1992, 1994; Covell, 1994. These methods are better ours in their generality because they require no knowledge of either actual or predicted secondary structure. On the other hand it’s difficult to see secondary structure at the resolution of these lattice models. In addition, often secondary structure is fitted onto the conformations at the end of construction. In computational performance, for small proteins these approaches compare very favorably with ours - a simplified representation which assigns a lattice point to every second residue can exhaustively sample shapes of proteins of up to 100 residues (Hinds & Levitt, 1992). However, because of an exponential increase in the number of shapes of a self-avoiding walk on a tetrahedral lattice, it’s difficult at this point to see a generalization of the lattice methods that would apply to larger molecules without a significant sacrifice in resolution.

One final method we want to mention is that of assembling structures from fragments of existing folds used by Simons et al., 1997, 1999. We find this method very appealing. It uses the information in the sequence to bias the assembly much the same way we use the ’patch’ information to bias helix-helix orientation. It is reported to produce best structures as good as the lowest RMSDs in our set. In addition, this method - as do the lattice models and the loop angle search - works on beta and mixed alpha/beta proteins: something we cannot yet do. Because we have not had the opportunity to test this method ourselves, we cannot comment on how its efficiency compares with ours. Our feeling is that sequence information in ’contact patches’ has more correlation on the global geometry of the structure than the local sequence/structure combination employed in these methods. Also our method is much faster than any of the other currently available sampling methods.

### C Advantages and disadvantages

In the beginning of his section we introduced two desirable qualities of a sampling method: speed, and the ability to produce structures close to the native fold. Our technique fills both of these requirements.

The speed of the method depends on the size of the protein and, more specifically, on the number of helices that we are trying to arrange. Each entry in table I takes from a fraction of a second to a few hours to produce. A naive count would estimate a factorial growth in the number of graphs, and an additional geometric growth in the number of structures sampled in each graph; however the actual increase is much less. For most structures of 4 helices or more, cutting the branches that violate self-avoidance results in a significant reduction of the number of conformations we need to sample. In our tests each additional helix increased the time of the runs by approximately a factor of 30. The set of conformations for our largest structure, 1fha, took approximately 4 hours to generate.

In addition to being fast, our method scales well - the use of branch-cutting helps because larger proteins are very constrained by self-intersection. Because of this our technique is able to sample molecules comparable in size to small domains, which are the largest single-chain structures one wishes to predict.

We chose to re-arrange helices using a ’patch’ contact database derived from existing structures. This approach significantly enriches the sampled ensemble with native-like structures. In addition, the graph-theoretical enumeration of relative orientations ensures that we sample all plausible regions of conformation space.

In its current incarnation our method also possesses some drawbacks. The main disadvantage of our technique is the requirement of a specified secondary structure. Ideally a search technique should also sample alternate secondary structure assignments. Currently we can only do this by specifying different assignments at the beginning of the procedure. Another significant drawback is the absence of a complimentary loop-building method. We need a method with the speed of a loop-library lookup methods, yet able to get a 3 Å or better approximation of the native loop (so that our near-native cores remain near-native). A final challenge to our approach is the lack of an obvious generalization to beta and alpha/beta proteins. The definitions of a sub-segment and the contact patch have to be significantly revised to adapt to beta sheets.

## V Future directions

The deficiencies of our technique, outlined in the previous paragraph, point the way to future developments. We are currently working on a fast loop building procedure for short (3-10 residues) loops. Our next project is to enhance sampling of alpha proteins to include variations of the boundaries of the helices.

A slightly more distant goal is the extension of the presented technique to construct $`\beta `$-sheet and mixed $`\alpha /\beta `$ proteins. In addition we are working on a fast preliminary discrimination function which would be used prior to the reconstruction of all atoms for each conformation.

## VI Acknowledgements

Many thanks to the present and former members of the Levitt group for help and discussions. We wish to thank Patrice Koehl for providing excellent routines for superimposing structures. Thanks to Prof. Joseph Rudnick for hints on graph theory. B.F. wishes to thank the A.P. Sloan Foundation, and the U.S. Department of Energy for financial support. This work was supported in part by grant DE-FG03-95ER62135 to M.L. from the U.S.Department of Energy.

Table I summarizes the sampling of helical cores for 25 proteins. Column 3 shows the number of residues that are assigned coordinates as well as the total length of the protein. Column 4 lists the range of CA RMSD of assigned coordinates from native coordinates. Columns 5 and 6 show the percentage of total structures that are closer than, respectively, 3$`\AA `$ and 6$`\AA `$ to the native structure. The last column displays the efficiency of our decoy generation method. It shows the ratio of the number of structures one needs to generate randomly in order to produce the best CA RMSD (computed using Eqn. 1) to the number of conformations in our ensemble.

Table II displays the results of pruning the decoy ensembles with disulfide bond information. The columns show shifts in the average Ca RMSD, proportions with RMSD below 3$`\AA `$ and 6$`\AA `$, and the reduction in the number of structures in each set.

Table III summarizes the results for decoy sets with all atoms’ coordinates assigned. The columns are similar to table I. The number of decoys in each ensemble is smaller than in table I because of pruning with a statistical pairwise potential. 2fha had no conformations below 6$`\AA `$, however 0.25 % of structures were closer than 10$`\AA `$ to the native.


text:
# Outflow from YSO and Angular Momentum Transfer

## 1. Introduction

Fragmentation of the clouds (cloud cores) is thought to be affected by the amount of angular momentum in it. On the angular momentum, there has been a problem that the specific angular momentum of new born stars (e.g. T Tauri stars) $`j_{}6\times 10^{16}(R_{}/2R_{})^2(P/10\mathrm{d}\mathrm{a}\mathrm{y})^1\mathrm{cm}^2\mathrm{s}^1`$ is much smaller than that of the parent molecular cloud cores $`j_{\mathrm{cl}}5\times 10^{21}(R/0.1\mathrm{pc})^2(\mathrm{\Omega }/4\mathrm{k}\mathrm{m}\mathrm{s}^1\mathrm{pc}^1)\mathrm{cm}^2\mathrm{s}^1`$. Here we assumed the observed velocity gradient comes from the rotational motion. Further, the specific orbital angular momentum of a binary whose separation is about $`R100`$AU is equal to $`j_{\mathrm{orbit}}4\times 10^{19}(R/100\mathrm{A}\mathrm{U})^{1/2}(M/1M_{})^{1/2}`$ $`\mathrm{cm}^2\mathrm{s}^1`$. If the angular momentum is conserved in the star forming process, neither binary nor single stars are formed.

Angular momentum is reduced by the effect of magnetic fields (B-fields). Considering an oblate spheroidal cloud is threaded by poloidal B-fields parallel to the symmetry axis, the rotation motion generates torsional Alfvén waves in the inlercloud medium, which propagate away from the cloud. The torsional Alfvén wave causes the intercloud medium to rotate. In this case, the angular momentum in the cloud decreases with a timescale as $`t_B\frac{\mathrm{\Omega }}{\dot{\mathrm{\Omega }}}\frac{\sigma }{2\rho _aV_A}=\frac{2\pi G\sigma }{B_0}(4\pi G\rho _a)^{1/2},`$ where $`\mathrm{\Omega }`$, $`\sigma `$, $`\rho _a`$, $`V_A`$, and $`B_0`$ represent, respectively, the angular rotation speed, the column density of the cloud, the ambient density, the Alfvén speed in the ambient medium, and the magnetic flux density outside the cloud (Ebert et al. 1960). Since the first factor $`2\pi G\sigma /B_0`$ is larger than unity for magnetically supercritical clouds, the timescale $`t_B`$ is comparable with the free-fall timescale in the low-density ambient medium and thus the angular momentum transfer may seem inefficient as long as the cloud contracts dynamically.

Molecular outflow begins at some epoch in the contraction stage of the cloud. The molecular outflow ejects matter and thus angular momentum. Consider the situation that 90 % of the inflowing matter forms a star and the rest 10 % is ejected by the molecular outflow. On the other hand, 0.1% of the angular momentum goes to the newly formed star and the rest is ejected by the outflow. If this occurs, the specific angular momentum of a newly formed star is reduced by a factor $`10^3`$ from that of the parent cloud. To confirm whether this mechanism works in a process of the star formation, we have studied the evolution of a slowly rotating cloud and the formation of molecular outflows.

## 2. Numerical Method & Model

We assume the gas is well coupled with the B-fields. The basic equations to solve are the ideal MHD equations with cylindrical symmetry and the Poisson equation for the gravitational potential. Since the gravitational collapse proceeds heavily nonhomologously, we have to cover the wide dynamic ranges for both the density and the size. This was attained by using the nested grid method, in which a number of grids with different spatial resolution are prepared (here we used 24 grids from the coarsest L0 to the finest L23) and calculated simultaneously. The finer grids are for small-scale structure near the core while the coarser ones cover the global structure.

As the initial condition, we assume a cylindrical isothermal cloud in hydrostatic balance. B-field strength is taken as proportional to $`\rho ^{1/2}`$ and the cloud’s rotation axis coincides with the axis of the cylinder. Density perturbations with $`\lambda `$ equal to that of the most unstable mode are added.

## 3. Run-away Collapse Stage

Isothermal clouds experience the run-away collapse, in which the central density increases very rapidly in a finite timescale. This is common for spherical clouds without rotation and B-fields (Larson 1969) and for a disk formed in a rotating isothermal cloud (Norman et al 1980). The gas falling down along the B-fields forms a contracting disk (pseudo-disk) even for rotating and magnetized clouds. Contraction continues if the equation of state remains isothermal. In this phase, there is no evidence for the outflow (Fig.1a,c).

## 4. Accretion Stage

After the central density reaches $`\rho _\mathrm{A}10^{10}\mathrm{cm}^3`$, the thermal radiation from dust grains is trapped in the central core and an adiabatic core is formed (Larson 1969). Since the dynamical timescale of the core becomes much shorter than that of the outer isothermal envelope, the isothermal gas begins to accrete to the nearly hydrostatic adiabatic atomic core (first core). To study the evolution of this accretion stage, we use a composite polytrope model as below $`\rho _\mathrm{A}`$ the gas is assumed isothermal; for $`\rho >\rho _\mathrm{A}`$ the gas becomes adiabatic and the polytropic index $`\mathrm{\Gamma }`$ is assumed as 5/3 (see Tomisaka 1998).

Figure 1 illustrates the change from the isothermal run-away collapse (Fig1.a, c) to the adiabatic accretion phase (Fig.1b, d). Figure 1a shows that B-fields (red lines) run almost perpendicularly to the disk (blue surfaces) in the run-away collapse phase. Toroidal B-fields and rotational motion are not predominant (Fig.1c). On the other hand, in the accretion phase B-fields are squeezed by the effect of inflow motion in the disk and they form an hourglass-shape configuration. Further, toroidal B-fields are formed by the rotational motions which are amplified in the accretion phase (Fig.1d). In this phase, outflow occurs (Fig.1b, d) at the distance $`50100`$ AU from the center. The time elapsed between these two is only about 1000 yr.

The Lorentz force $`𝐅=(c/4\pi )𝐣\times 𝐁`$ contains the toroidal component $`F_\varphi =(c/4\pi )𝐣_{\mathrm{pol}}\times 𝐁_{\mathrm{pol}}`$ when there exists the poloidal current ($`𝐣_{\mathrm{pol}}`$) which makes the toroidal B-fields $`B_\varphi `$. Therefore, it is shown that the magnetic torque $`N=rF_\varphi `$, which changes the distribution of the angular momentum, works effectively only after the toroidal B-field is amplified by the rotation motion. The toroidal B-fields become predominant after the core formation epoch. Further, Blandford & Peyne (1982) showed that efficient angular momentum transfer from the disk to the (cold) wind material occurs only when the angle between the B-field line and the disk is smaller than 60. This condition is satisfied only in the accretion phase (compare Fig.1a and.1b)

## 5. Angular Momentum Transfer

To study the effect of angular momentum transfer by the magnetic torque, specific angular momentum $`j`$ is plotted against the accumulated mass from the center $`M`$ in Figure 2 (Tomisaka 2000). The difference between the initial stage (dashed line) and the core formation epoch (dotted line) does not come from the angular momentum transfer but from the segregation of the mass with small $`j`$ to form the high-density central part. True transfer occurs after the core formation. The magnetic torque increases $`j`$ of the matter on the disk surface and the matter with excess $`j`$ is ejected as the outflow, while the magnetic torque decreases $`j`$ inside the disk and this leads to accelerate the accretion.

That is, B-fields changes the angular momentum distribution along one field line. As a result, in $``$7000 yr after

Figure 2: Specific angular momentum $`j`$ is plotted against the accumulated mass from the center $`M`$. Open circles, filled squares, and triangles denote, respectively, the beginning and the end of the run-away collapse phase and the accretion phase.

the core formation $`M0.2M_{}`$ is accumulated in the core and $`j`$ of the core becomes as small as $`10^{16}\mathrm{cm}^2\mathrm{s}^1`$ (solid line). Therefore, it is clearly shown that the angular momentum problem of newborn stars is solved taking into account the molecular outflow driven by the magnetic torque formed just outside the first core ($`R<100`$AU).

## References

Blandford, R.D., & Payne, D.G. 1982, MNRAS, 199, 883

Ebert, R., von Hoerner, S., & Temesvary, S. 1960, Die Entstehung von Sternen durch Kondensation diffuser Materie (Berlin: Springer), 315

Larson, R.B. 1969, MNRAS, 145, 271.

Norman, M.L., Wilson, J.R., & Barton, R.T. 1980, ApJ, 239, 968.

Tomisaka, K. 1998, ApJL 502, L163

Tomisaka, K. 2000, ApJL 528, L41; see also
http://yso.mtk.nao.ac.jp/$``$tomisaka/Presentation/iau200/main/index.htm