text:
# Neutrino mass patterns, WMAP, and Neutrinoless double beta decay

## I Introduction

Neutrino data collected in recent years almost established the hypothesis that neutrinos have masses and oscillate into each other, especially in the recent SNO sno  and KamLAND kamland  data. The SNO result sno  told us that the neutrino flux that arrived at the Earth consists of an electron neutrino component and a non-electron neutrino component, the sum of which matches in good agreement with the prediction of the solar standard model without oscillation ssm . The plausible explanation is that the electron neutrinos produced at the core of the sun oscillate into other (active) neutrino flavors, which are all detected in the SNO detector. On the other hand, KamLAND kamland  measured the anti-electron-neutrino fluxes collected from a number of nuclear plants. With these human-produced $`\overline{\nu }_e`$ fluxes at each plant and the initial flux known with a high precision, the solar neutrino problem was tested for the first time solely on Earth.

All the combined solar neutrino data (including KamLAND data) now favors the oscillation of the $`\nu _e`$ into a mixture of $`\nu _\mu `$ and $`\nu _\tau `$ in the large mixing angle (LMA) solution solar-fits

$$\mathrm{\Delta }m_{\mathrm{sol}}^27.1\times 10^5\mathrm{eV}^2;\mathrm{tan}^2\theta _{\mathrm{sol}}=0.45,$$
(1)
which gives $`\mathrm{sin}^22\theta _{\mathrm{sol}}0.86`$. The other solutions (small mixing angle, LOW, vacuum) are no longer favorable. In addition, the sterile neutrino is also strongly disfavored.

The atmospheric neutrino data from SuperK superk  and the data from K2K k2k  showed a maximal mixing between the $`\nu _\mu `$ and $`\nu _\tau `$ atm-fits  with a mass-square difference and a mixing angle:

$$\mathrm{\Delta }m_{\mathrm{atm}}^22.7\times 10^3\mathrm{eV}^2,\mathrm{sin}^22\theta _{\mathrm{atm}}=1.0.$$
(2)
Yet, another piece of important information came in very recently from WMAP wmap  on cosmic microwave background anisotropies. When combined with the data from 2dF Galaxy Redshift Survey, CBI, and ACBAR other , WMAP is able to pin down the amount of critical density contributed by the relativistic neutrinos. Thus, it puts an upper bound on the sum of neutrino masses weiler

$$\underset{i}{}m_i<0.71\mathrm{eV}$$
(3)
at 95% C.L. Such an upper limit already contradicts the result from LSND hitoshi .

There has been a claim of a positive signal observed in a neutrinoless double beta decay ($`0\nu \beta \beta `$) experiment 0nbeta . If it is real, it will add another constraint to the neutrino mass matrix. The positive signal in the $`0\nu \beta \beta `$ decay implies a nonzero $`m_{ee}`$ entry in the neutrino mass matrix arranged in the flavor basis:

$$\left(\begin{array}{ccc}m_{ee}& m_{e\mu }& m_{e\tau }\\ m_{e\mu }& m_{\mu \mu }& m_{\mu \tau }\\ m_{e\tau }& m_{\mu \tau }& m_{\tau \tau }\end{array}\right).$$

The 95% allowed range of $`m_{ee}`$ is $`0.110.56`$ eV with a best value of $`0.39`$ eV. 111 In Ref. 0nbeta , the allowed range of $`m_{ee}`$ could be widen to $`0.050.84`$ eV if allowing a more conservative estimate on the nuclear matrix elements. P. Vogel in PDG pdg  obtained a higher range for $`m_{ee}=0.41.3`$ eV using a different set of nuclear matrix elements. It also implied that the electron neutrino is majorana in nature. There are, however, arguments against this claim: see Ref. vissani ; aal , as well as counter arguments new .

Such a large value for $`m_{ee}\sqrt{\mathrm{\Delta }m_{\mathrm{atm},\mathrm{sol}}^2}`$, if true, gives a nontrivial modification to the neutrino mass patterns and mixings. A lot of possible mass textures that were proposed to explain the atmospheric and solar neutrino deficits become incompatible with the new $`0\nu \beta \beta `$ data. A well-known example is the Zee model zee  that can generate a bi-maximal mixing between the $`\nu _\mu `$ and $`\nu _\tau `$ and between the $`\nu _e`$ and the mixture of $`\nu _\mu `$ and $`\nu _\tau `$ jarl ; paul1 ; otto , however the Zee model guarantees the diagonal mass matrix elements to be zero. Therefore, it cannot explain the nonzero $`m_{ee}`$ implied by the $`0\nu \beta \beta `$ signal.

More and more neutrino data are accumulated since a few years ago. They all contribute to the understanding of neutrino masses and mixings. They are very important clues to the underlying theory for neutrino masses, and perhaps to other fermions as well. It is therefore very important to pin down the necessary patterns for neutrino mass matrix that can explain all the observations. The purpose of this short note is to derive the neutrino mass patterns from the most updated solar and atmospheric neutrino fits, and the WMAP constraint in Eq. (3). Finally, we show that it is rather marginally consistent with the claim in the $`0\nu \beta \beta `$ experiment.

Neutrino mass patterns under various theoretical or phenomenological assumptions or anatz have been considered in a lot of work xg .

## II Neutrino mass Patterns

To obtain the mass pattern it is easy to work bottom-up. We work in the basis where the charged-lepton mass matrix is diagonal such that the mixing information is entirely contained in the neutrino mass matrix. From the data we can write down specific forms of the mixing matrix $`U`$ and the diagonal mass matrix $`M_D`$ in the mass basis, then we can obtain the mass matrix $`M_\nu `$ in the flavor basis as

$$M_\nu =UM_DU^T,$$
(4)
where the flavor basis and mass basis are related by

$$\left(\begin{array}{c}\nu _e\\ \nu _\mu \\ \nu _\tau \end{array}\right)=\left(\begin{array}{ccc}U_{e1}& U_{e2}& U_{e3}\\ U_{\mu 1}& U_{\mu 2}& U_{\mu 3}\\ U_{\tau 1}& U_{\tau 2}& U_{\tau 3}\end{array}\right)\left(\begin{array}{c}\nu _1\\ \nu _2\\ \nu _3\end{array}\right).$$
(5)
Assuming no CP violation a general form of $`U`$ is given by

$$U=U_{23}U_{13}U_{12},$$
(6)
where $`U_{ij}`$ is the rotation matrix about the $`i`$ and $`j`$’th mass eigenstates. By considering the various data sets the general form of $`U`$ can be reduced to a simple form. Since we know that the atmospheric neutrino requires a maximal mixing between the $`\nu _\mu `$ and $`\nu _\tau `$, $`U_{23}`$ is given by

$$U_{23}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ 0& \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right).$$

Since the angle $`\theta _{13}`$ is constrained by CHOOZ chooz  to be small, we simply set $`U_{13}=I`$. For the rotation between the $`(1,2)`$ states we used a generic $`U_{12}`$ as

$$U_{12}=\left(\begin{array}{ccc}c& s& 0\\ s& c& 0\\ 0& 0& 1\end{array}\right),$$

where $`c=\mathrm{cos}\theta _{12},s=\mathrm{sin}\theta _{12}`$. Therefore, $`U`$ is given by

$$U=\left(\begin{array}{ccc}c& s& 0\\ \frac{s}{\sqrt{2}}& \frac{c}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ \frac{s}{\sqrt{2}}& \frac{c}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right).$$
(7)
For the mass matrix in the mass basis there are two cases: (i) normal and (ii) inverted mass hierarchies, with a mass scale $`m_0`$ ($`m_0`$ is considerably larger than the atmospheric mass scale) for the lightest one. In the normal hierarchy, the solar oscillation is between the two lighter states, while the atmospheric oscillation is between the heaviest and the lighter states. Recall the convention that the solar oscillation is between the “1” and “2” states, we put the diagonal mass matrix as

$$M_D^{\mathrm{normal}}=\left(\begin{array}{ccc}m_0& 0& 0\\ 0& (m_0+\delta )e^{i\varphi }& 0\\ 0& 0& (m_0+m)e^{i\varphi ^{}}\end{array}\right),$$
(8)
where $`\delta \mathrm{\Delta }m_{\mathrm{sol}}^2/(2m_0)`$ , $`m\mathrm{\Delta }m_{\mathrm{atm}}^2/(2m_0)`$, and we take the majorana phases $`\varphi `$ and $`\varphi ^{}`$ as 0 for the moment. We shall later show the results for other values of $`\varphi `$ and $`\varphi ^{}`$. Using Eqs. (4), (7), and (8), we obtain the neutrino mass matrix in the flavor basis:

$$M_\nu ^{\mathrm{normal}}=\left(\begin{array}{ccc}m_0+\delta s^2\hfill & \frac{cs}{\sqrt{2}}\delta \hfill & \frac{cs}{\sqrt{2}}\delta \hfill \\ \frac{cs}{\sqrt{2}}\delta \hfill & m_0+\frac{m}{2}+\frac{c^2}{2}\delta \hfill & \frac{m}{2}+\frac{c^2}{2}\delta \hfill \\ \frac{cs}{\sqrt{2}}\delta \hfill & \frac{m}{2}+\frac{c^2}{2}\delta \hfill & m_0+\frac{m}{2}+\frac{c^2}{2}\delta \hfill \end{array}\right).$$
(9)
For the case of inverted mass hierarchy the solar neutrino oscillation is between the two heavier states while the atmospheric neutrino oscillation is between the lightest and the heavier states. Recall the convention again that the solar oscillation is always between the “1” and “2” states, we used a $`M_D^{\mathrm{inverted}}`$:

$$M_D^{\mathrm{inverted}}=\left(\begin{array}{ccc}m_0+m+\delta & 0& 0\\ 0& m_0+m& 0\\ 0& 0& m_0\end{array}\right),$$
(10)
in which we have chosen the phase angles to be zero. Using Eqs. (4), (7), and (10) we obtain

$$M_\nu ^{\mathrm{inverted}}=\left(\begin{array}{ccc}m_0+m+\delta c^2\hfill & \frac{cs}{\sqrt{2}}\delta \hfill & \frac{cs}{\sqrt{2}}\delta \hfill \\ \frac{cs}{\sqrt{2}}\delta \hfill & m_0+\frac{m}{2}+\frac{s^2}{2}\delta \hfill & \frac{m}{2}+\frac{s^2}{2}\delta \hfill \\ \frac{cs}{\sqrt{2}}\delta \hfill & \frac{m}{2}+\frac{s^2}{2}\delta \hfill & m_0+\frac{m}{2}+\frac{s^2}{2}\delta \hfill \end{array}\right).$$
(11)
In general, we can use the data to fit the values of $`m_0,m,\delta ,\theta _{12}`$:

$`m_0m_{ee};m\mathrm{\Delta }m_{\mathrm{atm}}^2/(2m_0);`$
$`\delta \mathrm{\Delta }m_{\mathrm{sol}}^2/(2m_0),\theta _{12}=\theta _{\mathrm{sol}}.`$ (12)
We have already used a particular form of $`U`$ such that a maximal mixing is always between the “2” and “3” states.

Using the WMAP constraint we obtain

$`{\displaystyle \underset{i}{}}m_i`$ $``$ $`3m_0<0.71\mathrm{eV}`$ (13)
$``$ $`m_{ee}=m_0<0.24\mathrm{eV}.`$
Immediately, we can see that the $`m_{ee}`$ so obtained can only be in the lower mass range of the $`0\nu \beta \beta `$ decay data. In fact, the best fit value of 0.39 eV in the $`0\nu \beta \beta `$ decay is inconsistent with the WMAP result.

### II.1 Non-zero Majorana phases

Next we are going to show the results when we allow nonzero phases for $`\varphi `$ and $`\varphi ^{}`$. For illustrations we use (ii) $`\varphi =\pi ,\varphi ^{}=0`$, (iii) $`\varphi =0,\varphi ^{}=\pi `$, and (iv) $`\varphi =\pi ,\varphi ^{}=\pi `$. The case (i) $`\varphi =\varphi ^{}=0`$ has already been shown above. The neutrino mass matrices in the flavor basis for normal and inverted mass hierarchies are given, respectively, by

$`(\mathrm{ii})`$ $`\varphi =\pi ,\varphi ^{}=0:`$ (21)
$`\left(\begin{array}{ccc}\mathrm{cos}2\theta _{12}m_0\delta s^2& \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta & \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0+\frac{m_0+m}{2}\frac{c^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0\frac{m_0+m}{2}\frac{c^2}{2}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0\frac{m_0+m}{2}\frac{c^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0+\frac{m_0+m}{2}\frac{c^2}{2}\delta \end{array}\right)`$
$`\left(\begin{array}{ccc}\mathrm{cos}2\theta _{12}(m_0+m)+\delta c^2& \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta & \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)+\frac{m_0}{2}+\frac{s^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)\frac{m_0}{2}+\frac{s^2}{2}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)\frac{m_0}{2}+\frac{s^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)+\frac{m_0}{2}+\frac{s^2}{2}\delta \end{array}\right)`$
$`(\mathrm{iii})`$ $`\varphi =0,\varphi ^{}=\pi :`$ (29)
$`\left(\begin{array}{ccc}m_0+\delta s^2& \frac{cs}{\sqrt{2}}\delta & \frac{cs}{\sqrt{2}}\delta \\ \frac{cs}{\sqrt{2}}\delta & \frac{m}{2}+\frac{c^2}{2}\delta & m_0+\frac{m}{2}+\frac{c^2}{2}\delta \\ \frac{cs}{\sqrt{2}}\delta & m_0+\frac{m}{2}+\frac{c^2}{2}\delta & \frac{m}{2}+\frac{c^2}{2}\delta \end{array}\right)`$
$`\left(\begin{array}{ccc}m_0+m+\delta c^2& \frac{cs}{\sqrt{2}}\delta & \frac{cs}{\sqrt{2}}\delta \\ \frac{cs}{\sqrt{2}}\delta & \frac{m}{2}+\frac{s^2}{2}\delta & m_0+\frac{m}{2}+\frac{s^2}{2}\delta \\ \frac{cs}{\sqrt{2}}\delta & m_0+\frac{m}{2}+\frac{s^2}{2}\delta & \frac{m}{2}+\frac{s^2}{2}\delta \end{array}\right)`$
$`(\mathrm{iv})`$ $`\varphi =\pi ,\varphi ^{}=\pi :`$ (37)
$`\left(\begin{array}{ccc}\mathrm{cos}2\theta _{12}m_0\delta s^2& \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta & \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0\frac{m_0+m}{2}\frac{c^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0+\frac{m_0+m}{2}\frac{c^2}{2}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}m_0\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0+\frac{m_0+m}{2}\frac{c^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}m_0\frac{m_0+m}{2}\frac{c^2}{2}\delta \end{array}\right)`$
$`\left(\begin{array}{ccc}\mathrm{cos}2\theta _{12}(m_0+m)+\delta c^2& \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta & \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)\frac{m_0}{2}+\frac{s^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)+\frac{m_0}{2}+\frac{s^2}{2}\delta \\ \frac{\mathrm{sin}2\theta _{12}}{\sqrt{2}}(m_0+m)\frac{cs}{\sqrt{2}}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)+\frac{m_0}{2}+\frac{s^2}{2}\delta & \frac{1}{2}\mathrm{cos}2\theta _{12}(m_0+m)\frac{m_0}{2}+\frac{s^2}{2}\delta \end{array}\right)`$
The present data cannot distinguish the above patterns because they all give correct $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$, $`\mathrm{sin}^22\theta _{\mathrm{atm}}`$, $`\mathrm{\Delta }m_{\mathrm{sol}}^2`$, $`\mathrm{tan}^2\theta _{\mathrm{sol}}`$, and $`m_{ee}`$ by adjusting the parameters $`m_0,m,\delta ,\theta _{12}`$. The preference of one pattern over the others is a matter of taste. We suggest a fine-tuning criteria: the stability of the maximal mixing in the atmospheric sector and the large mixing in the solar sector provided by the mass matrix. The mass matrices in (i) are not stable, which is explained as follows. Let us look at the $`(2,3)`$ block of the matrix, which is of the form $`\left(\begin{array}{cc}1& ϵ\\ ϵ& 1\end{array}\right)`$, where $`ϵ1`$. The maximal mixing between the “2” and “3” states is very unstable. Once $`m_{\mu \mu }`$ and $`m_{\tau \tau }`$ differ very slightly, the maximal mixing is destroyed. We have verified that if $`m_{\mu \mu }`$ or $`m_{\tau \tau }`$ increases by 2% the resulting $`\mathrm{sin}^22\theta _{\mathrm{atm}}0.6`$, already out of the allowed range. On the other hand, in case (iii) the $`(2,3)`$ block is of the form $`\left(\begin{array}{cc}ϵ& 1\\ 1& ϵ\end{array}\right)`$. This is a rather robust form such that even when $`m_{\mu \mu }`$ or $`m_{\tau \tau }`$ is changed appreciably, the maximal mixing remains. We have checked even if $`m_{\mu \mu }`$ or $`m_{\tau \tau }`$ increases by a factor of two the maximal mixing solution remains for the atmospheric neutrino. This is easy to understand. As long as $`m_{\mu \tau }m_{\mu \mu },m_{\tau \tau }`$ we have the maximal mixing in the $`(2,3)`$ block. Unfortunately, for a similar reason the matrix in case (iii) is not stable with respect to the solar sector. Therefore, only cases (ii) and (iv) remain relatively stable with respect to both the solar and atmospheric sectors with a small change in the matrix elements.

In cases (ii) and (iv), the effective mass involved in $`0\nu \beta \beta `$ decay is

$$m_{ee}m_0\mathrm{cos}2\theta _{12}0.38m_0<0.09\mathrm{eV}$$
(38)
In these favorable cases, the upper limit on $`m_{ee}`$ is more stringent than that in Eq. (13). For this value of $`m_{ee}`$ it is barely inside the most conservative 95% C.L. range of the $`0\nu \beta \beta `$ experiment (on the other hand, well outside the range obtained by P. Vogel pdg .)

For a generic majorana phase $`\varphi `$ the effective mass $`m_{ee}`$ is bounded by

$$m_{ee}m_0\sqrt{c^4+s^4+2c^2s^2\mathrm{cos}\varphi }<(0.24\mathrm{eV})\sqrt{c^4+s^4+2c^2s^2\mathrm{cos}\varphi }.$$
(39)
## III Conclusions

We have derived a general neutrino mass pattern based on the most updated fits to the atmospheric and solar neutrino data. Using the WMAP constraint the common mass scale $`m_0`$ must be less than 0.24 eV, which implies an upper limit on the effective neutrino mass in $`0\nu \beta \beta `$ decay: $`m_{ee}<(0.24\mathrm{eV})\sqrt{\mathrm{cos}^4\theta _{12}+\mathrm{sin}^4\theta _{12}+2\mathrm{cos}^2\theta _{12}\mathrm{sin}^2\theta _{12}\mathrm{cos}\varphi }`$, where $`\varphi `$ is one of the majorana phases. The largest allowed value for $`m_{ee}=0.24`$ eV, which is already in the lower side of the allowed range indicated by the $`0\nu \beta \beta `$ experiment, and the best value of the experiment is already inconsistent with the $`m_{ee}`$.

In our favorable cases by considering the stability of the neutrino mass matrix, the $`m_{ee}`$ is further limited to $`m_{ee}<0.09`$ eV, which is just barely inside the most conservative range indicated by the experiment. It is well below the range obtained by P. Vogel pdg  using a different set of nuclear matrix elements.

We conclude that the WMAP constraint has a very important impact on the neutrino mass pattern and present a challenge to the Heidelberg-Moscow experiment. It also implies a sensitivity requirement on the future neutrinoless double beta decay experiments.

Note added: Similar conclusions on $`0\nu \beta \beta `$ decay have been obtained in Refs. weiler ; hitoshi . Other work related to the WMAP result can be found in Ref. josh .

This research was supported in part by the National Center for Theoretical Science under a grant from the National Science Council of Taiwan R.O.C.


format: text

text:
# The Formation of Stellar Clusters in Turbulent Molecular Clouds: Effects of the Equation of State

## 1 INTRODUCTION

Stars form alone and in groups, associations, and clusters (see Pudritz 2002 and Ward-Thompson 2002 for recent reviews). Yet the origin of this diversity in the process of star formation remains poorly understood. Recent work suggests that stars form by gravitational collapse of shock-compressed density fluctuations generated by the supersonic turbulence generally observed in molecular clouds (e.g. Elmegreen 1993; Padoan & Nordlund 1999; Klessen, Heitsch, & Mac Low 2000; Ostriker, Stone, & Gammie 2001; see Mac Low & Klessen 2003 for a recent review). If so, then star formation is controlled by turbulent fragmentation and dynamical collapse of individual, Jeans-unstable, protostellar clumps.

Fragmentation of gas clouds has been studied for more than a century, yet the process remains poorly understood. There have been analytical models (Jeans 1902; Low & Lynden-Bell 1976), and numerical investigations of the effects of various physical processes on the collapse, such as geometry and rotation of the clouds (see Bonnell & Bastien 1993 for a review), and magnetic fields (Galli et al. 2001). Recently, the effects of turbulence have been studied extensively in a series of three-dimensional simulations using both grid-based, Eulerian and particle-based, Lagrangian hydrodynamics (Klessen, Burkert, & Bate 1998; Klessen & Burkert 2000, hereafter Paper I; Klessen, Heitsch, & Mac Low 2000, Paper II; Heitsch, Mac Low, & Klessen 2001, Paper III; and Klessen 2001, Paper IV). These papers show that turbulence strongly influences the fragmentation of molecular clouds. Fast, clustered collapse and star formation occurs in regions with turbulence insufficient to support against gravitational collapse, while slow, scattered star formation results from strong turbulent support. However, these results are based on isothermal models with polytropic exponent $`\gamma =1`$ only. The role of the equation of state (EOS), which is essential in understanding the physical structure and stability of the turbulent clouds, remained unexplored.

The balance of heating and cooling in a molecular cloud can be approximately described by a polytropic EOS

$$P=K\rho ^\gamma ,$$
(1)
where $`K`$ is a constant, and $`P,\rho `$ and $`\gamma `$ are thermal pressure, gas density and polytropic exponent, respectively. The density structure generated by supersonic motions in highly compressible turbulence depends on the EOS. The stiffness of the EOS can largely control the density probability distribution function (PDF) in strongly compressible turbulence (Scalo 1998; Passot & Vázquez-Semadeni 1998; Vázquez-Semadeni & Garcia 2001). In particular, one-dimensional simulations by Passot & Vázquez-Semadeni (1998) show the shape of the density PDF varying with $`\gamma `$. They found that the PDF is log-normal when $`\gamma =1`$ (isothermal), and develops a power-law tail for $`\gamma 1`$. Spaans & Silk (2000) showed that radiatively cooling gas can be described by a piecewise polytropic EOS, in which the polytropic exponent $`\gamma `$ changes with gas density $`\rho `$. They conducted a detailed analysis including chemistry, thermal balance and radiative transfer, and found that in the interstellar medium the effective polytropic exponent $`\gamma `$ has a range of $`0.2<\gamma <1.4`$. They further predicted that the initial mass function (IMF) of protostellar cores would vary with $`\gamma `$, yielding a peaked IMF for $`\gamma >1`$ and a power-law function for $`\gamma <1`$.

In this paper, we perform three-dimensional simulations of self-gravitating, supersonic turbulence over the range of $`\gamma `$ identified by Spaans & Silk (2000) to investigate the effect of varying the EOS on the dynamical evolution and fragmentation behavior of turbulent clouds and the subsequent formation of protostellar clusters. To quantify the dependency on $`\gamma `$, we consider a true polytropic EOS, with $`\gamma `$ being strictly constant in each simulation regardless of the density, rather than the piecewise polytrope used by Spaans & Silk (2000) and earlier workers.

In §2 we describe our computational methods; in §3 we give results on fragmentation, the resulting mass distribution of protostellar cores, and the gas density PDF; and in §4 we conclude and speculate about the implications or our work for the formation of star clusters and primordial star formation.

## 2 METHODS

In analytical work on the stability of a self-gravitating, isothermal medium, Jeans (1902) derived a relation between the oscillation frequency $`\omega `$ and the wavenumber $`k`$ of small perturbations,

$$\omega ^2c_s^2k^2+4\pi G\rho =0$$
(2)
where $`c_s`$ is the isothermal sound speed, $`G`$ the gravitational constant, and $`\rho `$ is the gas density. The medium is unstable to fragmentation at all wavelengths greater than a critical length $`\lambda _J=2\pi /k_J`$, or equivalently all masses exceeding the Jeans mass,

$$M_J=\rho \lambda _J^3=\left(\frac{\pi }{G}\right)^{3/2}\rho ^{1/2}c_s^3$$
(3)
will collapse under their own weight. Note that we use a cubic definition of $`M_J`$.

In a polytropic cloud with EOS given by equation 1, the sound speed is

$$c_s=\left(\frac{dP}{d\rho }\right)^{1/2}=(K\gamma )^{1/2}\rho ^{(\gamma 1)/2},$$
(4)
so the Jeans mass will be given as a function of $`\gamma `$ and $`\rho `$,

$$M_J=\left(\frac{K\pi }{G}\right)^{3/2}\gamma ^{3/2}\rho ^{\frac{3}{2}(\gamma \frac{4}{3})}.$$
(5)
The specific internal energy is $`u=K\rho ^{\gamma 1}/(\gamma _{ad}1)`$, where $`\gamma _{ad}`$ is the adiabatic index, with $`\gamma _{ad}=5/3`$ for ideal gas. For isothermal gas, $`\gamma =1`$, and $`K=c_s^2`$.

We carry out simulations to determine when and how fragmentation occurs. We use a smoothed particle hydrodynamics (SPH) code (Benz 1990) in order to resolve several orders of magnitude in density during the collapse of shock-compressed regions. Our version of the code includes periodic boundary conditions (Klessen 1997) and can replace high-density cores with sink particles (Bate, Bonnell, & Price 1995). Sink particles accrete surrounding gas particles while conserving mass and momentum, but they only interact gravitationally. They prevent time step from becoming prohibitively short in very dense regions. This allows us to follow the dynamical evolution of the system over many free-fall times.

We set up a periodic cube with side length $`L=2`$ and total mass of unity such that initial density $`\rho _0=0.125`$, and sound speed $`c_s=0.1`$. The treatment of the polytropic EOS follows Bonnell & Bastien (1991), where we set $`K=(R_g/\mu )T_0\rho _0^{1\gamma }`$. We choose values consistent with the sound speed of $`R_g=2/3`$, $`\mu =1`$, and $`T_0=0.015`$ for the gas constant, mean molecular weight, and initial temperature of the cloud. Uniform turbulence is driven with the method described by Mac Low et al. (1998) and Mac Low (1999), adding energy over a small range of wavenumbers $`k`$. Driving in two wavenumber ranges is considered here, $`1k2`$ and $`7k8`$, corresponding to models $`B1h`$ and $`B3`$ in Paper II. In each case, the driving strength is chosen to ensure the same rms velocity despite the different driving wavenumbers, so that the average turbulent Jeans mass $`M_J_{turb}=3.2`$ and Jeans Number $`N_J=64`$, are the same. In a self-gravitating medium, the maximum resolvable density is determined by the requirement that the local thermal Jeans mass be resolved by the SPH kernel (Bate & Burkert 1997); our models all have the same resolution limit.

The models presented here are computed in normalized units. If scaled to mean densities of $`n(\mathrm{H}_2)10^2`$cm-3, and a temperature of $`11.4`$K (i.e. an isothermal sound speed $`c_s=0.2`$km s-1), values appropriate for a dark cloud like Taurus (Hartmann 1998), then our simulation cube holds a mass of $`M4\times 10^3`$M and has a size of $`L9`$pc. For details on the scaling relations, see Paper II.

We computed models with $`0.2\gamma 1.4`$, covering the range suggested for interstellar gas by Spaans & Silk (2000). We vary $`\gamma `$ in steps of 0.1 in otherwise identical simulations for both driving wavenumbers. Each simulation starts with a uniform density $`\rho _0=0.125`$. Driving begins immediately, while self-gravity is turned on at t=2.0, after the turbulence is fully established (see Paper II). Our models use $`\mathrm{200\hspace{0.17em}000}`$ particles. Twenty-six simulations were performed simultaneously using a serial version of the code on single 1 GHz Pentium III processors of the Parallel Computing Facility of the American Museum of Natural History (AMNH). Each computation took about 6 months.

## 3 RESULTS

### 3.1 Turbulent Fragmentation

Figure 1 and Figure 2 show the density distribution of the gas. To make the figures, a 2563 grid was filled with densities properly computed from the SPH kernel, and the result is displayed in three-dimensional projection (top row) and in a slice through maximum density one free-fall time after self-gravity is turned on for selected $`\gamma `$ and both driving wavenumbers. As found in previous studies (e.g. Mac Low 1999), the driving wavenumber strongly influences the density distribution. Driving with $`k=1`$–2 produces strong filamentary structure, while the density distribution for $`k=7`$–8 remains more uniform at large scale.

We find that the value of $`\gamma `$ is as important as the driving in determining collapse behavior. As $`\gamma `$ increases, the number of collapsed cores replaced by sink particles decreases, and the cores cluster less. For the same $`\gamma `$, driving with $`k=1`$–2 produces both far more, and more clustered dense cores than driving with $`k=7`$–8. The dense cores tend to collapse predominantly in filaments or at intersections of filaments. At $`\gamma 1.1`$, no fragmentation occurs by the time shown for driving with $`k=7`$–8.

The isothermal case, $`\gamma =1.0`$, agrees well with previous results in Papers II and III. Our current investigation agrees with the results in Paper I–III, that collapse tends to form clusters with high efficiency in regions with weak turbulence, while in regions with strong turbulence, sparse, slow collapse occurs. We now add the additional criterion that collapse and fragmentation depends strongly on $`\gamma `$.

In regions with $`\gamma <1`$, fragmentation occurs earlier and more frequently, while in regions where $`\gamma >1`$, fragmentation is retarded and less frequent.

Figure 3 compares the number of collapsed cores at different $`\gamma `$ for models with different driving wavenumbers. The rate at which new protostellar cores form differs for different $`\gamma `$. Models with low $`\gamma `$ form cores quickly, while models with large $`\gamma `$ form new cores more rarely. Again we see more cores for driving with $`k=1`$–2 than for $`k=7`$–8.

Rouleau & Bastien (1990) found an empirical scaling relation $`J_c\gamma ^2`$ of the critical Jeans number $`J_c=M/M_J`$ for collapse for elongated clouds. The Jeans number should be related to the number of fragments. However, with higher resolution and more statistics for a wide range of $`\gamma `$ in our simulations, we find something more like the inverse, but strongly dependent on the details of the driving. Fig 4 shows the relation between the value of $`\gamma `$ and the number of fragments at one free-fall time. The number of fragments drops quickly as $`\gamma `$ increases for driving with $`k=1`$–2, and more slowly for driving with $`k=7`$–8, which produces fewer fragments to begin with.

Figure 5 shows the accretion histories (the time evolution of the combined mass fraction of all protostellar cores) for all the individual $`\gamma `$ cases and for both driving wavelengths. Note that much less mass is accreted when driving with $`k=7`$–8. We focus on the $`k=1`$–2 model. The lower $`\gamma `$ is, the earlier fragmentation occurs and the larger the number of fragments that form. It also takes longer for the clouds with high $`\gamma `$ to fragment. The slope of the accretion curves are roughly the same, but the number of sinks is very different for different $`\gamma `$ at the same time, which suggests that the average core mass is different for different $`\gamma `$. A low-$`\gamma `$ environment produces a large number of collapsed cores of relatively low mass, while $`\gamma >1`$ results in the formation of considerably fewer, but more massive, cores. Since the overall mass growth rate of collapsing cores is closely related to the expected star formation rate, our result suggests that the overall star formation rate is higher in low-$`\gamma `$ clouds than in high-$`\gamma `$ clouds. It remains unclear what the final amount of mass in stars will be, as higher mass stars may drive surrounding gas away more quickly (see also Vázquez-Semadeni, Ballesteros-Paredes, & Klessen 2003).

We do not find any signs of collapse in the model with $`\gamma =1.4`$ even after several free-fall times. This can be explained with a Jeans mass analysis giving the critical condition for gravitational collapse. From Equation , we have

$$\frac{M_J}{\rho }\frac{3\gamma 4}{2}\rho ^{(3\gamma 6)/2}.$$
(6)
This implies that for $`\gamma <4/3`$ the Jeans mass $`M_J`$ decreases as density increases during collapse; for $`\gamma =4/3`$, $`M_J`$ remains constant; while if $`\gamma >4/3`$, then $`M_J`$ increases with density. Thus, any collapse will choke itself off when $`\gamma >4/3`$, as found in our numerical results.

Fragmentation is a complex process depending on the local conditions in the cloud, so the above analysis for an isolated spherical perturbation is an approximation. Indeed, fragmentation already ceases for $`\gamma >1.1`$ in the model driven with $`k=78`$, suggesting the importance of the details of the turbulence in determining collapse behavior.

### 3.2 Mass Spectrum

Observations suggest that the mass distribution of gas clumps in molecular clouds follows a power law, $`dN/dM=M^\nu `$, with typical values of the exponent being $`\nu 1.5`$ (e.g. Williams, Blitz, & McKee 2000). Salpeter (1955) derived a power law for the high mass stellar IMF in the same notation with $`\nu =2.35`$ (see also Scalo 1998; Kroupa 2002).

We measure clump-mass spectra in our models using a clump-finding algorithm similar to the one described by Williams, De Geus, & Blitz (1994) but working on all three spatial coordinates and adapted to make use of the SPH kernel smoothing procedure (for details see Appendix A in Paper I). Figure 6 shows the spectra of gas clumps and collapsed cores for models with different values of $`\gamma `$, driven with $`k=1`$–2. Three different evolutionary phases are shown: after turbulence has been fully established, but before self-gravity has been turned on, at $`t=2`$; and when the fraction of mass accumulated in collapsed cores (sink particles) has reached $`M_{}20\%`$ and $`40\%`$.

The mass distribution of both clumps and collapsed cores changes with $`\gamma `$, with the effect being most pronounced for the cores. In low-$`\gamma `$ models, the core mass spectrum at the high-mass end is roughly log-normal. As $`\gamma `$ increases, fewer but more massive cores emerge. When $`\gamma >1.0`$, the distribution is dominated by a few high mass cores, and the spectrum tends to flatten out. It is no longer fit by either a log-normal or a power-law. The clump mass spectra, on the other hand, do show power-law behavior on the high mass side, even for $`\gamma >1.0`$.

Our results suggest that massive stars can form in small groups or even in isolation in gas with $`\gamma >1.0`$. Spaans & Silk (2000) suggested that a stiff EOS with $`\gamma >1.0`$ should lead to a peaked IMF, biased toward massive stars, while an EOS with $`\gamma <1.0`$ results in a power-law IMF, in general agreement with our simulations.

### 3.3 Density Probability Distribution Function

Figure 7 shows mass and volume weighted PDFs of gas density $`p_m`$ and $`p_v`$ before self-gravity has been turned on for selected values of $`\gamma `$ and both driving wavenumbers. We also show Gaussian fits to the PDFs together with their error $`ϵ=|pp_{fit}|/p`$, summed over all points above $`10\%`$ of the peak values.

The mass-weighted density PDF $`p_m(\rho )`$ is calculated directly from the local density associated with each SPH particle, as described by Klessen (2000). The volume-weighted density PDF $`p_v(\rho )`$ is calculated on a cubic grid derived from the SPH density field by taking a cube of $`256^3`$ cells and using the SPH kernel smoothing procedure to find the density at the center of each cell. Both $`p_m(\rho )`$ and $`p_v(\rho )`$ are normalized so that $`_{\mathrm{}}^{\mathrm{}}p_m(\rho )𝑑\rho =_{\mathrm{}}^{\mathrm{}}p_v(\rho )𝑑\rho =1`$.

The density behind shocks in supersonic turbulence depends on the compressibility of the gas,

$$d\rho \frac{1}{\gamma }P^{\frac{1}{\gamma }1}dP,$$
(7)
So turbulent density fluctuations should increase as $`\gamma `$ decreases. Indeed, we see that in Figure 7, the width of $`p_m`$ and $`p_v`$ increases as $`\gamma `$ decreases.

In models driven with wavenumber $`k=78`$, there are super-Gaussian tails in the low-density end in the volume-weighted PDFs for all $`\gamma `$ and in the mass-weighted PDF for $`\gamma =0.2`$. One-dimensional simulations by Passot & Vázquez-Semadeni (1998) showed that the volume-weighted density PDF of supersonic turbulent gas displays a power-law tail at high densities for $`0<\gamma <1`$, becomes log-normal for $`\gamma =1`$ (isothermal), and develops a power-law tail at low density for $`\gamma >1`$, which is only in partial agreement with our results. The discrepancy may be due to the low resolution in low-density regions in our simulations, the difference between one- and three-dimensional simulations, or different driving.

A full investigation of the density PDFs under various conditions, with different codes, different dimensions (three-dimensional vs. one-dimensional), different resolution, different turbulence driving, different initial Mach number, and with and without magnetic field is beyond the scope of this paper; these issues will be addressed elsewhere.

#### 3.3.1 Comparison Between Mass and Volume Weighting

In order to compare the mass and volume weighted PDFs, we will find it useful to define $`s=\mathrm{ln}(\rho /\rho _0)`$; the mass and volume weighted PDFs of $`s`$ are then $`p_m(s)`$ and $`p_v(s)`$, respectively. We can relate $`p(s)`$ to $`p(\rho )`$ if we remember that for any montonic function $`y(x)`$, it can be shown that the PDF $`|p(y)dy|=|p(x)dx|`$. Therefore,

$$p_m(\rho )=p_m(s)\frac{ds}{d\rho }=\frac{p_m(s)}{\rho },$$
(8)
and

$$p_v(\rho )=p_v(s)\frac{ds}{d\rho }=\frac{p_v(s)}{\rho }.$$
(9)
Since $`p_m(\rho )dM/d\rho `$, and $`p_v(\rho )dV/d\rho `$ (Ostriker, Stone, & Gammie 2001), it follows that

$$p_m(\rho )\frac{dM}{dV}\frac{dV}{d\rho }\rho p_v(\rho ).$$
(10)
Then we can relate $`p_m(s)`$ to $`p_v(s)`$

$$p_m(s)=\rho p_m(\rho )=C\rho ^2p_v(\rho )=Ce^sp_v(s).$$
(11)
where C is a normalization constant.

If $`p_v(s)`$ follows a normal distribution,

$$p_v(s)=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{(ss_v)^2}{2\sigma ^2}\right),$$
(12)
where $`s_v`$ is the average (in volume) value of $`s`$, and $`\sigma `$ is the dispersion, then

$$p_m(s)=\frac{C}{\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{s_v^2(s_v+\sigma ^2)^2}{2\sigma ^2}\right)\mathrm{exp}\left(\frac{\left[s(s_v+\sigma ^2)\right]^2}{2\sigma ^2}\right).$$
(13)
Since $`_{\mathrm{}}^{\mathrm{}}p_m(s)𝑑s=1`$, the normalization

$$C=\mathrm{exp}\left(\frac{s_v^2(s_v+\sigma ^2)^2}{2\sigma ^2}\right).$$
(14)
Equation 13 can thus be rewritten as:

$$p_m(s)=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{(ss_m)^2}{2\sigma ^2}\right).$$
(15)
where $`s_m=s_v+\sigma ^2`$. This is also a normal distribution, with the same dispersion $`\sigma `$ as that of $`p_v`$ but shifted average value $`s_m`$. Using the volume PDF normalization $`\rho p_v(s)𝑑s=e^sp_v(s)𝑑s=1`$, it can be seen

$$s_m=s_v=\frac{\sigma ^2}{2}.$$
(16)
as noted, for example, by Ostriker, Stone, & Gammie (2001). This holds only if the density PDF is log-normal. (Note that in this case $`C=1`$.)

This derivation gives several properties of the density PDFs $`p_m`$ and $`p_v`$. If one of them is log-normal, so is the other; and if they are log-normal, the width $`\sigma `$ of both profiles should be equal (Eqs.  and ), and the peaks of the PDFs $`p_m(s)`$ and $`p_v(s)`$ should lie symmetrically around zero (Eq. ).

In Figure 7 the volume-weighted PDFs lie to the left of the mass-weighted PDFs, with $`s_ms_v`$, as predicted. The decent Gaussian fits to the peaks of the PDFs show that the assumption of log-normal behavior is not bad. Three-dimensional grid-based simulations of decaying turbulence by Ostriker, Stone, & Gammie (2001) show similar Gaussian PDFs.

A detailed examination of the moments of the density PDFs, shown in Fig 8, shows that these PDFs are, in fact, not perfect log-normals. The first moment (mean) shows that $`p_m`$ and $`p_v`$ of models with both driving wavenumbers $`k=12`$ and $`k=78`$ are not exactly symmetrically distributed. This could be due to turbulent intermittency in our small boxes producing deviations from a purely Gaussian distribution. The Gaussian behavior expected from a well sampled distribution is more visible in models driven with $`k=78`$ compared to $`k=12`$. In fact, $`p_m`$ and $`p_v`$ of $`k=78`$ are more equal in width, as seen from the second moment (variance). This is a sign of better sampling – more modes contribute to the overall velocity field, thus the central limit theorem is more appropriately applied to the $`k=78`$ than to the $`k=12`$ case (Klessen 2000).

#### 3.3.2 Effect of Self-Gravity

Figure 9 shows the evolution of the mass-weighted gas density PDF $`p_m`$ of the model driven with $`k=12`$ for different values of $`\gamma `$. The same three evolutionary phases as defined in §3.2 are shown: the initial distribution, and times at which dense cores have accreted $`M_{}20\%`$ and $`M_{}40\%`$ of the mass. These PDFs do not include gas accreted into the cores. We again show the best-fit Gaussian for the part of each PDF above 10% of peak.

During the dynamical evolution, the PDFs in the non-isothermal cases develop pronounced positive deviations from the the Gaussian fit, while the PDF in the isothermal case remains well fit by a Gaussian (of increasing width), consistent with the results in Klessen (2000). The high-density tails gradually diminish as $`\gamma `$ increases from 0.2 to 1.0, vanishing almost entirely at $`\gamma =1.0`$. At $`\gamma >1.0`$, though, the tails develop again and get stronger as $`\gamma `$ continues to increase. We do not fully understand why the tails only appear in non-isothermal cases.

The first four moments of the density PDFs for each of the three evolutionary phases discussed previously are shown in Figure 10. We can see that the moments vary as collapse proceeds, demonstrating that self-gravity plays an important role in shaping the density PDF.

The plot of the skewness shows that the PDFs of turbulence without self-gravity are always skewed to low densities, regardless of $`\gamma `$. Shock compression comes at a cost: the gas swept up in shock fronts has to come from somewhere, and those regions now contain only low-density gas. As the filling factor of shocks is low (Smith, Mac Low, & Heitsch 2000) most of the volume is occupied by these low-density regions resulting in PDFs with negative skewness. The situation changes, however, when self-gravity comes into play. Gravitational contraction produces high densities in addition to the turbulent shock compression, the PDF’s in the later stages of evolution are thus strongly biased towards high-density gas and the skewness becomes positive.

## 4 SUMMARY AND DISCUSSION

We have investigated the effect of varying the equation of state (EOS) on the formation of stellar clusters in turbulent molecular clouds using a set of three-dimensional numerical simulations. Our results show that the EOS plays an important role in the fragmentation of the clouds, the determination of the initial mass function (IMF) of the protostellar cores, and the shape of the gas density PDFs.

The ability of interstellar gas clouds to fragment under the action of self-gravity decreases with increasing polytropic exponent $`\gamma `$ in the range $`0.2<\gamma <1.4`$ relevant for Galactic molecular clouds. The total amount of mass in collapsed cores, however, appears to remain roughly constant through that range. At least half of this material is expected to accrete onto the protostellar system at the center of each collapsing core (e.g. Klessen 2001). Small values of $`\gamma `$ thus lead to the formation of dense clusters of predominantly low-mass stars, while $`\gamma >1`$ results in the formation of isolated and massive stars. Fragmentation and collapse ceases entirely for $`\gamma >1.4`$ as expected from analytic arguments.

The mass spectra of both clumps and collapsed cores change with $`\gamma `$, with that effect being most pronounced for the cores. For $`\gamma 1.0`$, their mass spectrum appears roughly log-normal. As $`\gamma `$ increases, fewer but more massive objects emerge. When $`\gamma >1.0`$, the distribution is dominated by high mass objects only, and the spectrum tends to flatten out. It is fit by neither a log-normal nor a power-law. The clump mass spectra, on the other hand, do show power-law behavior at the high mass end for all $`\gamma `$.

The density PDF changes with $`\gamma `$ as well. The width of the profile decreases as $`\gamma `$ increases, as expected from analytic arguments about shock compression. Both mass- and volume-weighted PDFs show imperfect log-normal distributions. Self-gravity helps shape the PDF. During dynamical collapse, the PDF develops a high-density tail, which is the imprint of local collapse of high density regions.

In a polytropic gas cloud, any temperature fluctuation corresponds to a density fluctuation, $`\mathrm{\Delta }\mathrm{log}_{10}(T)\mathrm{\Delta }\mathrm{log}_{10}(\rho )(1\gamma )0.6`$. We note that from the gas density PDFs in Figure 7, the maximum width of the profiles is about 12, which corresponds to $`\mathrm{\Delta }\mathrm{log}_{10}(\rho )5`$. So for the extreme case $`\gamma =0.2`$, $`\mathrm{\Delta }\mathrm{log}_{10}(T)3.5`$. For a cloud with an average temperature of about 10 K the maximum temperature of the low-density gas therefore could reach values up to $`10^4`$ K, if $`\gamma `$ is extremely low. This still falls in the temperature range observed in star forming regions.

The principal motivation of this paper is to investigate how the fragmentation behavior of a turbulent cloud and the subsequent formation of stars are affected by physical state of the cloud, as described by the EOS. Since $`\gamma `$ reflects the balance between heating and cooling, $`\gamma `$ is expected bo be low in regions of strong cooling, while $`\gamma `$ may be large where the cooling mechanism is insufficient (Spaans & Silk 2000). Unfortunately, the true chemical state and composition of interstellar gas is difficult to assess, so exact values of $`\gamma `$ are difficult to estimate observationally.

The possibility that isolated, massive stars form when $`\gamma >1`$ is of interest, not only because ionizing radiation, stellar winds, and supernovae from massive stars play an important role in the interstellar medium, but also because massive stars are rarely observed in isolation. Rather, they are usually found as members of rich stellar clusters. Recently, however, Lamers et al. (2002) has reported observations of massive stars found in isolation or associated only with very small groups of lower-mass stars in the bulge of M51. Massey (2002) also reported finding apparently isolated, massive, field stars in both the Large and Small Magellanic Clouds. Our simulations suggest that molecular cloud regions forming isolated, high-mass stars have $`\gamma >1`$. In this case turbulent cloud fragmentation tends to produce isolated, massive collapsing cores that will feed onto isolated high-mass stars. Yorke & Sonnhalter (2002) demonstrated that high-mass stars can form in collapsing gas clumps via disk accretion very much the same way as low mass stars are believed to do. There is thus no need for a dense cluster environment where high-mass stars could build up by collisions of lower-mass stars (as suggested by Bonnell, Bate, & Zinnecker 1998).

A theory of galactic-scale star formation partially based on the universality of log-normal density PDFs has been proposed by Elmegreen (2002). Our results do show a log-normal form at least for the peak of the PDFs, containing most of the gas. However, the dispersion of the PDF varies substantially as the EOS changes (Fig. 8), so that the fraction of gas available for star formation would depend significantly on the EOS in such a theory. This raises the question of how it could still be able to explain a global Schmidt Law across many different galaxies with different metallicities, radiation fields, interstellar pressures, and so forth.

High resolution simulations by Abel, Bryan, & Norman (2000, 2002) of the formation of the first star suggest that fragmentation during the collapse of metal-free pregalactic halos is rather inefficient, resulting in the formation of single, massive stars rather than clusters of lower-mass objects. In the absence of metals, inefficient cooling may result in high $`\gamma `$, perhaps helping to explain these results.

We do note, however, that simulations by Bromm, Coppi, & Larson (1999, 2002) show greater multiplicity. They use similar chemistry, but examine more massive, more isolated halos, with substantially larger limiting mass resolution. The EOS will have less effect on these larger scales, which are more dominated by gravity. If the metal-free gas indeed shows a high effective $`\gamma `$, our models would suggest that their collapsing regions will show little further fragmentation if followed down to stellar mass scales.

We thank M. Fall, H. Lamers, H. Zinnecker, and S. Glover for valuable discussions, and D. Janies and W. Wheeler for their work on the AMNH Parallel Computing Facility, which we used for the computations presented here. We are grateful to the anonymous referee for useful comments and corrections. YL thanks the AIP for its warm hospitality, and the Kade Foundation for support of her visits there. M-MML acknowledges partial support by NASA ATP grant NAG5-10103, and by NSF CAREER grant AST99-85392. RSK acknowledges support by the Emmy Noether Program of the Deutsche Forschungsgemeinschaft (DFG, KL1358/1).


text:
# Interferometric Detection of Linear Polarization from Sagittarius A* at 230 GHz

## 1 Introduction

Linear and circular polarization are important diagnostics of synchrotron sources. In active galactic nuclei and galactic microquasars, polarization is useful for characterizing particle energy, magnetic field configuration and jet composition (Marscher & Gear, 1985; Wardle et al., 1998; Fender et al., 2002). Variability in polarization is important for diagnosing shocks and the stability of accretion processes.

Here we present new observations of linear and circular polarization in Sagittarius A\*, the compact nonthermal radio source in the Galactic Center. Sgr A\* is associated with a $`3\times 10^6\mathrm{M}_{\mathrm{}}`$ black hole (Melia & Falcke, 2001). While much is known about this source, there is still dispute over the nature of accretion and outflow in the source. Advection dominated accretion flows (ADAFs), convection dominated accretion flows (CDAFs), Bondi-Hoyle accretion, Blandford-Königl-type jets and hybrids of these ideas have been used to characterize the spectrum and size of Sgr A\* (e.g., Narayan et al., 1998; Quataert & Gruzinov, 2000; Liu & Melia, 2001; Falcke & Markoff, 2000; Yuan et al., 2002). In general models account for the low luminosity of Sgr A\* through a very low accretion rate, through a very low radiative efficiency or through a combination of these. Estimates of the accretion rate vary by seven orders of magnitude. Linear polarization observations can break this degeneracy.

These new observations complement past work that has shown a complex picture of the polarization in Sgr A\*. The emission shows no linear polarization between 1.4 and 112 GHz at limits of 0.1 to 2% (Bower et al., 1999a, c, 2001b). Surprisingly, the emission is circularly polarized between 1.4 and 43 GHz (Bower et al., 1999b; Sault & Macquart, 1999; Bower et al., 2002a). The circular polarization is strong and variable, sometimes showing a highly inverted spectrum. Models account for the circular polarization and absence of linear polarization through conversion of linear to circular polarization and bandwidth depolarization of the linear polarization with cold electrons (Ruszkowski & Begelman, 2002; Beckert & Falcke, 2002). A model that involves plasma modes in the vicinity of a black hole has also been developed (Broderick & Blandford, 2000).

The polarization properties above 112 GHz are less certain. JCMT observations at 150, 220, 375 and 400 GHz (2.0, 1.35, 0.85 and 0.75 mm wavelength, respectively) indicate that Sgr A\* is polarized at a level of $`10\%`$ (Aitken et al., 2000, hereafter A00). However, these observations suffer from the poor angular resolution of the JCMT (22′′ at 220 GHz). A substantial amount of polarized and unpolarized radiation must be subtracted from the flux detected at the Galactic Center to estimate the polarization of Sgr A\*. The apparent sharp transition in linear polarization apparent between 112 and 150 GHz from these observations is very unusual and has important implications for the source physics (Quataert & Gruzinov, 2000; Agol, 2000; Melia et al., 2000).

We present interferometric observations of the linear polarization in Sgr A\* at 230 GHz (1.3 mm wavelength). Made with the BIMA array, these observations are approximately two orders of magnitude higher resolution than those of the JCMT. Thus, we are able to eliminate the effects of extended polarized and unpolarized emission and isolate the emission of Sgr A\* itself. We present the observations in §2, discuss the results in §3 and summarize in §4.

## 2 Observations and Data Reduction

Data were obtained on 4 epochs (Table 1) in the B and C configurations of the BIMA array (Welch et al., 1996). The measured single sideband system temperature on Sgr A\* ranged from 1000 to 2800 K for the B-configuration data, and 1500 to 9000 K for the C-configuration data. Both upper (230.7 GHz) and lower sideband (227.3 GHz) were sampled with a bandwidth of 800 MHz. The data were obtained from switched polarization observations (see Bower et al., 1999c). The polarization on each antenna is switched between left and right circular polarization (L and R) using a Walsh switching pattern of period 16. Using an integration time of 11 second, all polarization combinations (LL, RR, LR and RL) are sampled every 5 minutes.

The data were phase self-calibrated by using the LL and RR polarizations averaged over 5 minutes. Phase calibration was applied to all 4 polarizations LL, RR, LR, and RL. Amplitude calibration of the BIMA array is determined from measurements of the antenna gains ($`150`$ Jy K-1) using planet observations at a number of epochs. The antenna gains are stable and no further amplitude calibration was applied.

The observed amplitudes are sensitive to atmospheric coherence and to antenna pointing errors. Both errors affect polarized and unpolarized emission equally, implying that the polarization fractions are a more reliable measurement than the total or polarized intensities. The antenna pointing was checked on the quasar J1733-1304 and Uranus. We estimate the magnitude of these errors from measurements of the flux density for J1733-1304 in February through May 2002. The measurements are $`2.2\pm 0.2`$ Jy at 0.85 mm (JCMT); $`2.3\pm 0.3`$ (BIMA) and $`3.0\pm 0.3`$ (OVRO), at 1.3 mm; and, $`4.3\pm 0.3`$ (BIMA) and $`3.7\pm 0.3`$ (OVRO) at 3 mm (B. Matthews and J. Carpenter, private communication). These are consistent with a coherence loss $`<30\%`$ for our 1.3 mm polarization measurements.

Polarization calibration was applied using the polarization leakage terms derived from observations in the same mode of 3C 273 and 3C 279 on 02 February 2002. The polarization leakage for each antenna ranges from $`1`$ to 5% on all but one antenna. The polarization results for Sgr A\* were identical with and without antenna 8, which had leakage terms $`10\%`$. Comparison of the polarization calibrations in 8 epochs over the past 1.5 years shows that the leakage terms are stable and reproducible to $`1.5\%`$ rms (Bower et al., 2002b). This leads to a random systematic error in the polarization of $`0.5\%`$.

We compared the polarized flux density for 3C 279 measured with BIMA at 86 and 230 GHz and measured with the VLA at 22 and 43 GHz over the past 2 years (Taylor & Myers, 2000; Bower et al., 2002b). The fractional polarization shows the same behavior from 22 to 230 GHz. The VLA and BIMA position angles differ by $`10`$ degrees. It is not known whether this difference is due to source structure or an intrinsic error in the VLA or BIMA calibration.

The polarization leakage corrections were applied by averaging the switched polarization data over intervals ranging from 2 to 30 minutes. We calibrated upper and lower sideband observations separately for each epoch. The errors in the polarization were consistent with the thermal noise for averaging times from 5 to 30 minutes.

Combining all the data weighted by system temperature and spatially filtered to include only (u,v)-data between 20 and 300 $`k\lambda `$, we obtain the images shown in Figure Interferometric Detection of Linear Polarization from Sagittarius A\* at 230 GHz. With robust weighting to minimize the noise from sidelobes and thermal noise we obtain a synthesized beamwidth 3.6 x 0.9 arcsec with a theoretical rms noise level 10 mJy. We combined the data in several combinations to look for changes in the polarization with time, frequency and array configuration. The results are summarized in Table 2. Fluxes were found by fitting point sources at the center of the field. LSB and USB refer to the lower and upper sideband data sets including all epochs. B, C1, C2 and C3 refer to the four epochs summarized in Table 1. C refers to a combination of all three C-configuration observations. For the combined result, we find the linear polarization to be $`7.2\pm 0.6\%`$ in position angle $`139\pm 4`$ degrees.

## 3 Discussion

### 3.1 Dust Polarization

These observations confirm the existence of high fractional linear polarization in or near Sgr A\* at 230 GHz. They can be used to eliminate the possibility of dust contamination of the polarized signal. A00 estimated that their central $`22^{\prime \prime }`$ beam included 3.55 Jy of (unpolarized) free-free emission and 0.75 Jy of dust emission at 220 GHz which left $`2.2\pm 0.5`$ Jy for Sgr A\*. A00 assume that the dust emission is 3% polarized at a position angle of 100 degrees based on the larger scale properties of the dust polarization. If we assume that the dust and free-free emission are distributed smoothly over this beam, then we expect a factor of $`\frac{22\times 22}{3.6\times 0.9}=150`$ decrease in these flux densities for the BIMA result. Thus, even 100% polarized dust emission could produce only $`5`$ mJy of the observed 89 mJy of polarized flux density.

Could the dust be significantly clumped to produce the observed polarization? If all of the dust emission were concentrated in our central pixel, then the minimum dust polarization fraction is $`12\%`$. However, this leaves only 0.5 Jy flux density for Sgr A\* itself, which is very low given the known spectrum and history of millimeter flux density observations. The minimum flux density observed for Sgr A\* at 1.3 mm is $`>1`$ Jy (Zhao et al., 2002). Introducing a decorrelation factor $`\xi 1`$ gives a flux density of $`0.5\xi ^1`$ Jy for Sgr A\* and a dust polarization fraction $`12\xi ^1\%`$. Our limit on gain errors derived from J1733-1304 places the constraint $`\xi >0.7`$, which doesn’t substantially increase the flux density for Sgr A\*. Interferometric observations of dust polarization indicate that dust polarization rarely exceeds 20% and only in cases of dense star-forming molecular clouds (Rao et al., 1998; Lai et al., 2002). An unlikely coincidence would be necessary to place such a rare, highly clumped and highly polarized cloud at the location of Sgr A\*.

We conclude therefore that the effects of dust polarization on our measurement are negligible and that the emission is associated with Sgr A\* itself.

### 3.2 Variations with Wavelength and Time

While the polarization fraction we measured agrees with that of A00, the position angle of $`139\pm 4`$ degrees does not. A00 measured a total polarization of 4.1% (266 mJy) in position angle $`89\pm 3`$ degrees in the central pixel at 220 GHz. Some of this polarized emission is dust which A00 assume to be polarized at 3% (25 mJy) in position angle 100 degrees. We can estimate the actual dust parameters necessary to make our measurements agree. Since our measured Sgr A\* polarization is almost orthogonal in $`(Q,U)`$-space to the central pixel of A00, the required dust polarization for agreement is extreme: 300 mJy (or 40% of the estimated 750 mJy of dust flux density) in position angle 80 degrees. This high polarization fraction is inconsistent with other measured polarization fractions as discussed above. Thus, we cannot easily reconcile the A00 220 GHz result with ours through the dust properties.

We can estimate the rotation measure from the position angle data as a function of frequency (Figure Interferometric Detection of Linear Polarization from Sagittarius A\* at 230 GHz). The lack of bandwidth depolarization in the BIMA results places a strong upper limit of $`10^8\mathrm{rad}\mathrm{m}^2`$ on the RM. The LSB and USB results for BIMA are nearly identical, placing an upper limit of $`2\times 10^6\mathrm{rad}\mathrm{m}^2`$. The best fit to these data points is $`3\times 10^5\mathrm{rad}\mathrm{m}^2`$ but the error is $`>10^6\mathrm{rad}\mathrm{m}^2`$. Although we can reconcile all the BIMA and JCMT position angles with RM $`3\times 10^6\mathrm{rad}\mathrm{m}^2`$ and three phase wraps between 150 and 400 GHz, a RM this large would lead to substantial depolarization in the broad (40 GHz) JCMT bandpasses. We would expect to find a much higher polarization fraction in the BIMA results which have much narrower bandpasses. In fact, based on the nearly equal polarization fractions in the BIMA and JCMT results, we can conclude again that the polarization must be $`<10^6\mathrm{rad}\mathrm{m}^2`$.

However, if we exclude the A00 220 GHz result, we find that all of the remaining polarization detections can be explained as the result of a single RM without leading to substantial depolarization. We measure $`\mathrm{RM}=4.3\pm 0.1\times 10^5\mathrm{rad}\mathrm{m}^2`$ and a zero-wavelength position angle of $`181\pm 2`$ degrees. Including the A00 220 GHz result gives $`4.3\times 1.6\times 10^5\mathrm{rad}\mathrm{m}^2`$. This error could be construed as including the effects of systematic measurement error and variability. If we exclude both the 150 and 220 GHz A00 result, we find $`\mathrm{RM}=4.4\pm 0.4\times 10^5\mathrm{rad}\mathrm{m}^2`$. While the exclusion of the low frequency JCMT points in the RM fit appears ad hoc, note that these are the least reliable of the JCMT results due to the much larger beam size. Nevertheless, this RM measurement must be considered tentative.

We cannot eliminate variability as an explanation for the differences in position angle that are seen. The A00 results were obtained in three separate epochs spread over 5 months nearly three years earlier than the BIMA results. However, we detect no change in the position angle of the linear polarization with time in the BIMA data (Figure Interferometric Detection of Linear Polarization from Sagittarius A\* at 230 GHz). The reduced $`\chi ^2`$ is 0.5 for the four epochs under the hypothesis of a constant position angle. Further, the A00 measurements at 375 and 400 GHz were separated by 5 months but show no apparent change. On the other hand, the fractional linear polarization in the BIMA data is slightly variable. The final epoch C3 is significantly different from the mean. The reduced $`\chi ^2`$ is 2.0 for the four epochs under the hypothesis of a constant polarization.

We summarize detections and upper limits for the linear polarization fraction in Figure Interferometric Detection of Linear Polarization from Sagittarius A\* at 230 GHz. There is a sharp transition in the polarization fraction between 100 and 200 GHz. Our previous BIMA measurement that led to an upper limit at 112 GHz cannot be explained by bandwidth depolarization. Bandwidth depolarization occurs at 112 GHz in 800 MHz of bandwidth when $`\mathrm{RM}>1.0\times 10^7\mathrm{rad}\mathrm{m}^2`$.

On the other hand, angular depolarization, which requires that the RM change the position angle by $`>180`$ degrees across the face of the source, is marginally adequate. At 112 GHz, an RM a few times $`10^5\mathrm{rad}\mathrm{m}^2`$ is just sufficient. Thus, a fully turbulent accretion region with an outer scale comparable to the source size ($`<1\mathrm{A}\mathrm{U}`$, see discussion below) could depolarize the source. However, measured RM flucutuations in the Galactic Center and Cygnus regions have a maximum of $`<10^3\mathrm{rad}\mathrm{m}^2`$ on degree scales and are $`<<1\mathrm{rad}\mathrm{m}^2`$ when extrapolated to scales below 1 arcsec (Lazio et al., 1990; Yusef-Zadeh et al., 1997). Nevertheless, turbulence in the accretion environment may be substantially different from that in the interstellar medium.

An alternative explanation is that Sgr A\* consists of an unpolarized source and a polarized source with a steeply inverted spectrum that dominates the spectrum above 230 GHz. A two component model is natural for both inflow and outflow models since they must account for the sharply increasing spectrum at millimeter wavelengths (Serabyn et al., 1997; Falcke et al., 1998). All models (see next section) that fit the centimeter to submillimeter spectrum include at least two distinct populations of radiating particles. If the spectrum rises as $`\nu ^{2.5}`$, then the source would depolarize by a factor of 0.16 between 230 and 112 GHz, which is sufficient to account for the absence of detection at 112 GHz. Spectral indices of $`<1`$ and $`>4`$ are excluded by the results.

### 3.3 Circular Polarization

Circular polarization is not detected in any of these data sets. For the cumulative data set the circular polarization is less than the thermal rms noise of $`1\%`$. Only in the B array data is there a hint of a detection. Here the circular polarization is two times the noise.

### 3.4 Origin and Propagation of the Linear Polarization

The maximum observed RM in the Galactic Center is on the order of a few times $`10^3\mathrm{rad}\mathrm{m}^2`$ (Yusef-Zadeh et al., 1997). However, this occurs in a nonthermal filament approximately 0.5 degrees away from Sgr A\* and may not probe the densest ionized regions towards the Galactic Center. A hyperstrong interstellar scattering screen significantly broadens the image of Sgr A\* at all wavelengths (Lo et al., 1998; Lazio & Cordes, 1998; Bower et al., 2001a). The physical parameters of this scattering screen are not fully constrained. However, we can estimate an upper limit to the RM of the screen. The maximum RM is $`10^4\mathrm{rad}\mathrm{m}^2`$ for models of the ionized surfaces of molecular clouds (Yusef-Zadeh et al., 1994; Lazio & Cordes, 1998; Bower et al., 1999a). However, the RM may arise in a region different from the scattering region. Sgr A\* is embedded in or behind Sgr A West as well as a dense, ionized halo that extends over $`5^{}`$ (Pedlar et al., 1989; Anantharamaiah et al., 1999). This halo has density of $`10^210^3\mathrm{cm}^3`$ and a scale of 10 pc. For magnetic field strengths of a milliGauss, the maximum RM is $`10^7\mathrm{rad}\mathrm{m}^2`$. Similarly, Sgr A West has a density $`10^2\mathrm{cm}^3`$, a scale of 1 pc and milliGauss field strength, giving a maximum RM $`10^5\mathrm{rad}\mathrm{m}^2`$. Thus, while the interstellar scattering screen does not contribute significantly to the measured RM, material which closer to Sgr A\* but not directly associated with accretion and outflow could.

Millimeter VLBI observations require a size for Sgr A\* of 72 gravitational radii or less at 43 GHz and sizes decreasing as $`\nu ^1`$ at higher frequencies (Bower & Backer, 1998; Lo et al., 1998; Krichbaum et al., 1998; Doeleman et al., 2001). Theoretical temperature profiles indicate that the emission comes from a volume within a few gravitational radii of Sgr A\* (Melia & Falcke, 2001). In ADAF or CDAF models, the 230 GHz source is embedded within a quasi-spherical accretion region. Imaging by Chandra at X-ray wavelengths indicate that Sgr A\* is extended on a scale consistent with the radius of this larger region, suggesting a concentration of hot thermal plasma (Baganoff et al., 2001b).

The emission must propagate through the accretion region without substantial depolarization. Using the measured RM $`4\times 10^5\mathrm{rad}\mathrm{m}^2`$, we can estimate the accretion rate. (Note that the following conclusions also hold for the upper limit RM of $`2\times 10^6\mathrm{rad}\mathrm{m}^2`$ determined from the BIMA data). The implied accretion rates differ between models because of the different temperature versus radius relations. We can calculate the Bondi-Hoyle accretion rate assuming equipartition and a uniform magnetic field equation (5) from Bower et al. (1999a). To produce the measured RM, the accretion rate must be $`10^7\mathrm{M}_{\mathrm{}}\mathrm{y}^1`$. However, this accretion rate is too low to produce the observed spectrum by two orders of magnitude for both the ADAF and Bondi-Hoyle models. An accretion rate of $`10^5\mathrm{M}_{\mathrm{}}\mathrm{y}^1`$ for these models can model the spectrum but requires an RM $`10^{10}\mathrm{rad}\mathrm{m}^2`$. On the other hand, CDAF models can model the spectrum with $`\dot{M}10^8\mathrm{M}_{\mathrm{}}\mathrm{y}^1`$, which produces $`\mathrm{RM}10^6\mathrm{rad}\mathrm{m}^2`$ (Quataert & Gruzinov, 2000). Jet models require an accretion rate of $`10^8\mathrm{M}_{\mathrm{}}\mathrm{y}^1`$, which is consistent with a low RM (Falcke & Markoff, 2000). Melia et al. (2000) predict $`10^{11}\mathrm{M}_{\mathrm{}}\mathrm{y}^1`$, which significantly underpredicts the RM. This model requires that the RM originates external to Sgr A\*.

Two effects may mitigate the conclusions discussed above. These are magnetic field reversals and deviations from equipartition between particle and magnetic field energy. If the magnetic field reverses frequently, then the rms field in the region is a better indicator of the characteristic field for a RM calculation (Ruszkowski & Begelman, 2002). A model predicting a high RM for a uniform field may produce an RM lower by orders of magnitude if there are many field reversals. However, the number of field reversals necessary to account for the difference between the ADAF or Bondi-Hoyle models and our measured RM is on the order of $`10^8`$, which is extreme. Reducing the magnetic field strength substantially below equipartition values will reduce the RM by the same factor. However, these reductions will be offset by the necessary increase in particle density to model the same spectrum.

In jet models, the emission originates in a region of similar scale that can be located an arbitrary distance from the black hole. Jet emission is likely to originate in a compact shocked component near the base of the jet (Yuan et al., 2002). The compressed magnetic field will be perpendicular to the axis of the jet, implying that the jet is oriented North-South. Increasing magnetic field order towards the base of the jet can account for the sharply increasing polarization fraction with frequency. High frequency VLBI observations of quasars and BL Lac objects indicate that electric vector position angles are typically but not always aligned with the jet axis (Lister, 2001). The inclination angle of a jet in Sgr A\* is not well-constrained. However, an angle close to the line of sight is possible given that highly inclined sources often produce the highest linear polarization in compact radio sources. This raises the possibility that the radio/millimeter emission is relativistically beamed, which may also account for the absence of a visible jet or ejected components despite strong variability (e.g., Bower & Backer, 1998).

The jet model does not require propagation through the accretion medium. However, the high value for the RM suggests that it may pass through some of the accretion or outflow medium (Beckert & Falcke, 2002). In AGN cores, the maximum RMs measured are on the order of 2000 $`\mathrm{rad}\mathrm{m}^2`$ (Taylor, 2000). This argues that we are probing a region in Sgr A\* quite different from that in other AGN.

The RM we measure is not inconsistent with a lower limit for the RM from spectropolarimetry at 4.8 and 8.4 GHz (Bower et al., 1999a). No linear polarization was detected for RMs as high as $`10^7\mathrm{rad}\mathrm{m}^2`$ in this experiment. Different source physics may account for the absence of polarization at the longer wavelengths. This will result if the lower frequency source is intrinsically depolarized due to greater internal field disorder, for example. Alternatively, foreground beam depolarization across the much larger lower frequency source may account for the absence of linear polarization.

### 3.5 Predictions for Future Observations

An important expectation of emission from a jet is variable linear polarization (Marscher & Gear, 1985). Variable polarization may also be a measure of the stability of accretion in a disk model. For example, Liu & Melia (2002) have suggested that precession of the disk on a timescale of 100 days is a consequence of a very low value for the spin of the black hole. Such precession could lead to a change in the position angle of the polarization if the disk orientation is linked to the magnetic field structure. We have not seen convincing evidence for this position angle variation yet, except perhaps in the discrepancy between A00 and our result. Significant variability in the total intensity at centimeter, millimeter, submillimeter and X-ray wavelengths on timescales of hours to years has been reported (Wright & Backer, 1993; Young et al., 2001; Baganoff et al., 2001a; Zhao et al., 2001; Bower et al., 2002a). Correlating polarization variability with these events could be an important clue to the origin of the polarization. However, the stable polarization that we have measured over 60 days and that A00 measured over 5 months between 375 and 400 GHz suggests that such connections are possibly subtle or require intensive monitoring. If changes in the polarization are associated with X-ray flares, then they may have a timescale of hours.

The detection of a RM and a power law index for the polarization fraction can be tested with observations at other wavelengths. BIMA observations can be spread more widely over frequency space. The 1.3 mm receivers are sensitive from 210 to 270 GHz. With properly calibrated quarter wave plates, BIMA could measure the expected difference in position angle of $`20`$ degrees over this range. BIMA can also measure differences in the power law index for the polarization fraction over this same range. Interferometric 150 GHz observations should see a substantially decreased polarization fraction but a significantly different polarization angle than observed at 230 GHz. Higher frequency observations may see the polarization fraction continue to increase.

Variability in the RM can be used as a diagnostic of the accretion or outflow medium. Sgr A\* shows variability on all timescales from hours to years at radio wavelengths (Zhao et al., 2001; Bower et al., 2002a). This variability may be indicative of changes in the accretion or outflow rates or conditions, which may be reflected in the RM. Refractive scintillation timescales are on the order of 1 month at 230 GHz for a velocity of 100 $`\mathrm{kms}^1`$.

The detection of linear polarization raises the possibility of polarimetric VLBI. Images made at 230 GHz or higher frequencies will be only weakly affected by scattering and may show the effects of general relativity (Falcke et al., 2000; Broderick & Blandford, 2000; Bromley et al., 2001). This may be the best probe of the close environment of a black hole available.

## 4 Conclusions

We have confirmed the existence of strong linear polarization at millimeter wavelengths in Sgr A\* with high resolution observations by the BIMA array. These eliminate the possibility of confusion from polarized dust emission. The polarized emission probably comes from within 10 gravitational radii of the black hole. This is one of the closest probes of the environment of a black hole that exists.

Propagation of the polarized radiation through the accretion and outflow media provides important limits on the physical characteristics in these regions. The measured RM appears to favor a very low accretion rate which excludes standard ADAF models and Bondi-Hoyle models with very high accretion rates for Sgr A\*.

The existence of linear polarization at these high frequencies opens up several new methods for investigation of this source. Variations in polarization angle and fraction as well as rotation measure may provide critical constraints on models. Very long baseline interferometric imaging of linearly polarized emission in the vicinity of the black hole may reveal a wealth of effects, including those of general relativity.

The BIMA array is operated by the Berkeley-Illinois-Maryland Association under funding from the National Science Foundation.