text:
# Linear Polarization in Gamma-Ray Bursts: The Case for an Ordered Magnetic Field

## 1. Introduction

The first detection of polarization in an optical afterglow was in GRB 990510, where a degree of polarization $`P=1.7\pm 0.2\%`$ ($`1.6\pm 0.2\%`$) was measured at $`t_{\mathrm{obs}}18`$hr ($`21`$hr) after the burst.111See Covino et al. 2003a for references to the above observations as well as to subsequent polarization measurements. Since then, $`P13\%`$ has been detected in a few additional afterglows,222There is one exception, $`P=9.9\pm 1.3\%`$, measured in GRB 020405 at $`t_{\mathrm{obs}}=1.3`$days (Bersier et al. 2003). Significantly lower values ($`P1.52\%`$) were measured in this afterglow at $`t_{\mathrm{obs}}=1.2,2.2`$, and $`3.3`$days by other groups, with a similar $`\theta _p`$. If real, this behavior has no simple explanation in any of the existing models. some of which showed a temporal variation in $`P`$ (Rol et al. 2000; Barth et al. 2003), but typically the position angle (PA) $`\theta _p`$ showed little or no change. A few other afterglows produced only upper limits, $`P25\%`$.

The polarization is attributed to synchrotron emission behind a shock wave. It thus depends on the local magnetic field configuration, which determines the polarization at each point of the afterglow image, and on the global geometry of the shock, which determines how the polarization is averaged over the (unresolved) image. We make a distinction between magnetic field configurations that are axially symmetric about the normal $`\widehat{n}_{\mathrm{sh}}`$ to the shock surface and those that are not.333In this picture, the field is tangled over very small scales and possesses this symmetry when averaged over regions of angular size $`1/\gamma `$, where $`\gamma `$ is the Lorentz factor of the shocked fluid. The first category gives no net polarization for a spherical flow and is assumed in most previous works (Sari 1999; Ghisellini & Lazzati 1999; Medvedev & Loeb 1999; Granot et al. 2002; Rossi et al. 2002). These models take the field to be completely random in the plane of the shock, and the polarization is usually attributed to a jet viewed somewhat off axis. For a structured jet, $`P`$ has one peak near the jet break time $`t_j`$ (when $`1/\gamma `$ increases to $`\theta _0`$, the initial jet opening half-angle), whereas for a uniform jet $`P`$ has 2 (or 3) peaks near $`t_j`$, with $`P`$ passing through zero and $`\theta _p`$ changing by $`90^{}`$ between the peaks.

The second category can produce net polarization even for a spherical flow. One example is the “patchy coherent field” model of Gruzinov & Waxman (1999), where the observed region consists of $`N50`$ mutually incoherent patches of angular size $`1/\gamma `$, within each of which the field is fully ordered.444A similar model was used to study the linear polarization induced by microlensing of GRB afterglows (Loeb & Perna 1998). This model predicts $`PP_{\mathrm{max}}/N^{1/2}10\%`$ (where $`P_{\mathrm{max}}6070\%`$ is the maximum $`P`$ of local synchrotron emission in a uniform magnetic field) and simultaneous (random) variability in $`P`$ and $`\theta _p`$ on timescales $`\mathrm{\Delta }t_{\mathrm{obs}}t_{\mathrm{obs}}`$ . The idea behind this model is that $`1/\gamma `$ is the angular size of causally connected regions, and for a magnetic field that is generated in the shock itself this is the largest scale over which the field can be coherent.

However, if the magnetic field were ordered on an angular scale $`\theta _B1/\gamma `$, then the resulting $`P`$ could approach $`P_{\mathrm{max}}`$. Such a situation can be realized if an ordered field exists in the medium into which the shock propagates. For a typical interstellar medium, the postshock field would still be very weak (with the magnetic energy a fraction $`ϵ_B10^{10}`$ of the internal energy), but it would be higher ($`ϵ_B10^4`$; Biermann & Cassinelli 1993) if the shock expands into a magnetized wind of a progenitor star, and higher yet ($`ϵ_B0.010.1`$) if it propagates into a pulsar-wind bubble (PWB), as expected in the supranova model (Königl & Granot 2002). A strong ordered field component is likely to exist in the original ejecta and could give rise to a high value of $`P`$ ($`P_{\mathrm{max}}`$) in both the prompt GRB and the reverse-shock emission.

We calculate the polarization for a jet with a tangled magnetic field in § 2, investigate the effects of adding an ordered field component in § 3, and discuss the results in § 4.

## 2. Polarization from a Jet with a Tangled Magnetic Field

Synchrotron emission is generally partially linearly polarized. In terms of the Stokes parameters: $`V0`$, $`\theta _p=\frac{1}{2}\mathrm{arctan}(U/Q)`$, and $`P=(Q^2+U^2)^{1/2}/I`$. As the Stokes parameters are additive for incoherent emission, they can be calculated by summing over all the contributions from different fluid elements to the flux at a given observed time $`t_{\mathrm{obs}}`$. In practice, the flux is calculated by dividing $`t_{\mathrm{obs}}`$ into bins of size $`\delta t_{\mathrm{obs}}`$, and assigning to the appropriate time bin the contribution to $`F_\nu `$ from the emission of each 4-volume fluid element $`d^4x`$,

$$dF_\nu (t_{\mathrm{obs}},\widehat{n},𝐫,t)=\frac{(1+z)^2}{d_L^2}\frac{j_\nu ^{}^{}(𝐫,t)\mathrm{\Theta }(\delta t_{\mathrm{obs}}2|t\frac{\widehat{n}𝐫}{c}t_{\mathrm{obs}}|)d^4x}{\gamma ^2(𝐫,t)\left[1\widehat{n}𝐯(𝐫,t)/c\right]^2\delta t_{\mathrm{obs}}},$$
(1)
$$\left\{\begin{array}{c}U/I\\ Q/I\end{array}\right\}=\left(dF_\nu \right)^1dF_\nu \left\{\begin{array}{c}P\mathrm{sin}2\theta _p\\ P\mathrm{cos}2\theta _p\end{array}\right\},$$
(2)
where $`j_\nu ^{}^{}`$ is the local emissivity, $`𝐯`$ is the fluid velocity, $`\widehat{n}`$ is the direction to the observer, $`\mathrm{\Theta }`$ is the Heaviside step function and the summation in Eq. 2 is over $`𝐫`$ & $`t`$, for fixed $`t_{\mathrm{obs}}`$ & $`\widehat{n}`$.

In this section, we consider only a random magnetic field $`B_{\mathrm{rnd}}`$ that is tangled on angular scales $`1/\gamma `$, with axial symmetry w.r.t. $`\widehat{n}_{\mathrm{sh}}`$. The field anisotropy is parameterized by $`b2B_{}^2/B_{}^2`$, where $`B_{}`$ ($`B_{}`$) is the magnetic field component parallel (perpendicular) to $`\widehat{n}_{\mathrm{sh}}`$. The local polarization of the emission from a given fluid element is555For simplicity, we assume $`j_\nu ^{}^{}(B^{}\mathrm{sin}\chi ^{})^ϵ`$, where $`\mathrm{cos}\chi ^{}=\widehat{n}^{}\widehat{B}^{}`$, with $`ϵ=2`$ (Sari 1999), although in reality $`ϵ`$ may be different.

$$\frac{P_{\mathrm{rnd}}(\theta ^{})}{P_{\mathrm{max}}}=\frac{\left(B_{}^2B_{}^2/2\right)\mathrm{sin}^2\theta ^{}}{B_{}^2\mathrm{sin}^2\theta ^{}+(1+\mathrm{cos}^2\theta ^{})B_{}^2/2}=\frac{(b1)\mathrm{sin}^2\theta ^{}}{2+(b1)\mathrm{sin}^2\theta ^{}},$$
(3)
(Gruzinov 1999; Sari 1999), where $`\mathrm{cos}\theta ^{}=\frac{\mu v/c}{1\mu v/c}`$, $`\mu \widehat{n}\widehat{n}_{\mathrm{sh}}`$. For $`P>0`$, the polarization is in the direction of $`\widehat{n}\times \widehat{n}_{\mathrm{sh}}`$.

The magnetic field configuration behind the shock, and hence the value of $`b`$, cannot be easily deduced from first principles. It was suggested that small-scale postshock fields can be generated by a two-stream instability (e.g., Medvedev & Loeb 1999), which predicts $`b1`$. However, it is not clear whether the magnetic fields produced in this way survive in the bulk of the postshock flow (Gruzinov 1999, although see Frederiksen et al. 2003), or whether this is the dominant tangling mechanism. The relatively low observed values of $`P`$ ($`34\%`$) suggest that $`0.5b2`$ if the polarization is due to a jet with a shock-generated field.

Turbulence in the postshock region (possibly induced by a microinstability) could amplify and isotropize the field, keeping $`b`$ close to $`1`$. As each fluid element moves downstream from the shock transition, it is sheared by the flow. For a Blandford-McKee (1976) self-similar blastwave solution with an ambient density $`\rho _{\mathrm{ext}}r^k`$, the length of a fluid element in the directions parallel and perpendicular to $`\widehat{n}_{\mathrm{sh}}`$ scales with the self-similarity variable $`\chi `$ (where $`\chi =1`$ at the shock front and increases with distance behind the shock) as $`L_{}\chi ^{\frac{92k}{2(4k)}}`$ and $`L_{}\chi ^{\frac{1}{4k}}`$. Therefore, the stretching of each fluid element in the radial direction is larger than in the tangential direction. This would increase $`b`$ while maintaining axial symmetry about $`\widehat{n}_{\mathrm{sh}}`$. The relevant value of $`b`$ here is the average over the postshock region, weighted by the emissivity. If the turbulence only persists over a small distance ($`0<\chi 11`$, but still $``$ the plasma skin depth) then, since most of the emission is from $`\chi `$ a few, we may have $`1<b`$a few.666On the other hand, if an intrinsically isotropic turbulence persists over a large portion of the emission region, this could reduce the effect of shearing on $`b`$, resulting in $`0<(b1)1`$.

Figure 1 shows the polarization lightcurves for $`b=0`$ and $`\mathrm{}`$, based on the jet model of Kumar & Panaitescu (2000). For all viewing angles $`\theta _{\mathrm{obs}}\theta _0`$ there are two peaks in $`P`$, a little before and after $`t_j`$: $`P`$ passes through zero in between these peaks as $`\theta _p`$ changes by $`90^{}`$. This result is similar to that of Ghisellini & Lazzati (1999), who did not consider lateral spreading of the jet, and differs from that of Sari (1999), who assumed $`\theta _{\mathrm{jet}}(t_{\mathrm{obs}}`$$`>`$$`t_j)=1/\gamma `$, since the lateral spreading in the jet model that we use is smaller than the one used by Sari. The main distinction between the $`b<1`$ and $`b>1`$ cases is a $`90^{}`$ difference in $`\theta _p`$, but this prediction can only be tested if one can independently determine the direction from the l.o.s. to the jet axis. This may in principle be done by measuring the direction of motion of the flux centroid $`\widehat{n}_c`$ (Sari 1999): for $`b<1`$, $`\widehat{P}`$ is perpendicular to (aligned with) $`\widehat{n}_c`$ before (after) $`t_j`$, whereas for $`b>1`$ the situation is reversed. More generally, one would then be able to test if indeed $`\widehat{P}\widehat{n}_c`$ or $`\widehat{P}\widehat{n}_c`$, as expected for a pure $`B_{\mathrm{rnd}}`$ field, or if the angle between them is different, as would generally be the case if an ordered field component were also present.

## 3. The Effects of an Ordered Magnetic Field

We now add an ordered magnetic field component $`B_{\mathrm{ord}}`$ to the random field $`B_{\mathrm{rnd}}`$ considered in § 2. In passing through the shock transition, the parallel component of the ambient magnetic field $`B_{\mathrm{ext}}`$ remains unchanged but the transverse component is amplified by a factor equal to the fluid compression ratio, which for $`\gamma 1`$ is $`4\gamma `$. Thus typically $`B_{}B_{}`$ behind the shock. For simplicity we assume that $`B_{\mathrm{ord}}`$ lies in the plane of the shock and is fully ordered and that $`B_{\mathrm{ext}}`$ is uniform, so that $`B_{\mathrm{ord}}`$ is coherent over the entire shock.

It is most convenient to sum over the Stokes parameters associated with $`B_{\mathrm{rnd}}`$ and $`B_{\mathrm{ord}}`$ separately and combine them at the end.777This is valid in the limit where the two components are associated with distinct fluid elements. Alternative schemes for combining $`B_{\mathrm{rnd}}`$ and $`B_{\mathrm{ord}}`$ may produce a somewhat different polarization. The direction of polarization $`\widehat{P}_{\mathrm{ord}}`$ of the emission from the ordered component is perpendicular to its projection ($`𝐁_{\mathrm{ord}}`$) on the plane of the sky; $`\widehat{P}_{\mathrm{rnd}}`$ is either along the plane containing the l.o.s. and the jet symmetry axis (for $`P_{\mathrm{rnd}}>0`$) or perpendicular to that direction (for $`P_{\mathrm{rnd}}<0`$). The total polarization and the PA are given by

$`P`$ $`=`$ $`\left({\displaystyle \frac{\eta P_{\mathrm{ord}}}{1+\eta }}\right)\left[1+\left({\displaystyle \frac{P_{\mathrm{rnd}}}{\eta P_{\mathrm{ord}}}}\right)^22\left({\displaystyle \frac{P_{\mathrm{rnd}}}{\eta P_{\mathrm{ord}}}}\right)\mathrm{cos}2\delta \right]^{1/2},`$ (4)
$`\theta _p`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{arctan}\left({\displaystyle \frac{\mathrm{sin}2\delta }{\mathrm{cos}2\delta \eta P_{\mathrm{ord}}/P_{\mathrm{rnd}}}}\right),`$ (5)
where $`\eta I_{\mathrm{ord}}/I_{\mathrm{rnd}}B_{\mathrm{ord}}^2/B_{\mathrm{rnd}}^2`$ is the ratio of the observed intensities in the two components and $`\theta _p,\delta `$ are measured as illustrated in the upper left inset of Figure 1.

For $`B_{\mathrm{ord}}`$ we find that the PA as a function of the polar angle $`\theta `$ from the l.o.s. and the azimuthal angle $`\varphi `$ (measured from $`\widehat{B}_{\mathrm{ord}}`$) is given, in the relativistic ($`\gamma 1`$) limit, by $`\theta _p=\varphi +\mathrm{arctan}(\frac{1y}{1+y}\mathrm{cot}\varphi )`$, where $`y(\gamma \theta )^2`$. We have $`I_\nu =I_\nu ^{}^{}(\nu /\nu ^{})^3`$, with $`I_\nu ^{}^{}\nu ^\alpha [1(\widehat{n}^{}\widehat{B}_{\mathrm{ord}}^{})^2]^{ϵ/2}`$, where $`\nu /\nu ^{}\frac{2\gamma }{1+y}`$ and $`1(\widehat{n}^{}\widehat{B}_{\mathrm{ord}}^{})^2\left(\frac{1y}{1+y}\right)^2\mathrm{cos}^2\varphi +\mathrm{sin}^2\varphi `$. The Stokes parameters are given by $`(U,Q)/IP_{\mathrm{max}}=𝑑\mathrm{\Omega }I_\nu (\mathrm{sin}2\theta _p,\mathrm{cos}2\theta _p)/𝑑\mathrm{\Omega }I_\nu `$. For a spherical flow or a jet at $`t<t_j`$, when the edge of the jet is not visible, $`𝑑\mathrm{\Omega }=_0^{2\pi }𝑑\varphi _0^{y_{\mathrm{max}}}𝑑y`$, $`U=0`$ and $`\theta _p=\pi /2`$. For $`ϵ=2`$ and $`y_{\mathrm{max}}=1`$ we obtain $`Q/IP_{\mathrm{max}}f(\alpha )=\frac{(2^{4+\alpha }1)(2+\alpha )(3+\alpha )}{2^{7+\alpha }28+2\alpha [2^{3+\alpha }57+(2^{3+\alpha }1)\alpha ]}`$, so $`P_{\mathrm{ord}}/P_{\mathrm{max}}=f(\alpha )0.900.93`$ for $`0<\alpha <1.5`$. For $`\alpha =\frac{p1}{2}`$, $`P_{\mathrm{max}}=\frac{p+1}{p+7/3}`$, $`ϵ=\frac{p+1}{2}`$ (i.e. PLS G in Granot & Sari 2002), our analytic result corresponds to $`p=3`$, for which $`P_{\mathrm{ord}}/P_{\mathrm{max}}=\frac{93}{101}0.92`$ and $`P_{\mathrm{ord}}=\frac{279}{404}0.69`$. For $`y_{\mathrm{max}}1`$, $`ϵ=2`$ we obtain $`f(\alpha )=\frac{(2+\alpha )(3+\alpha )}{8+5\alpha +\alpha ^2}`$; for $`\alpha =\frac{p1}{2}`$ and $`p=3`$, $`P_{\mathrm{ord}}/P_{\mathrm{max}}=\frac{6}{7}0.86`$ and $`P=\frac{9}{14}0.64`$. The diference between the $`y_{\mathrm{max}}=1`$ and $`y_{\mathrm{max}}1`$ results may be relevant to the prompt GRB (see § 4), where the tail of a pulse (which corresponds to $`y>1`$) is predicted to be less polarized than its peak ($`y1`$). When the edge of the jet is visible, the limits of integration over $`d\mathrm{\Omega }`$ change. As this causes relatively small modifications in $`\theta _p`$ and $`P_{\mathrm{ord}}`$, we use the analytic expressions above for simplicity.

Figure 2 depicts a sample of polarization lightcurves in which both $`\eta `$ and $`b`$ are taken to be independent of time. In this case the $`B_{\mathrm{ord}}`$-induced polarization is constant (in both $`P`$ and $`\theta _p`$) throughout the afterglow. Interestingly, a similar polarization signature could be produced by dust in our galaxy or in the GRB host galaxy. In the latter case, however, the polarization would likely be accompanied by absorption that would redden the spectrum: this could in principle make it possible to estimate the level of the galactic dust contribution and thereby determine the fraction of such a constant-polarization component that is intrinsic to the source. Since $`P_{\mathrm{ord}}P_{\mathrm{max}}`$ whereas $`P_{\mathrm{rnd}}`$ is typically much smaller, we find that for $`\eta =1`$, and even for $`\eta =1/3`$, the polarized intensity is still dominated by $`B_{\mathrm{ord}}`$ ($`\eta P_{\mathrm{ord}}>P_{\mathrm{rnd}}`$), with $`B_{\mathrm{rnd}}`$ only inducing relatively small fluctuations around the $`B_{\mathrm{ord}}`$-induced values of $`P`$ and $`\theta _p`$. For $`\delta =45^{}`$ the fluctuations in both $`P`$ and $`\theta _p`$ are very small in this parameter range.

If $`B_{\mathrm{ord}}`$ dominates the polarization, then, by equation (4), the time evolution of $`P`$ follows that of the intensity ratio $`\eta `$. The low measured values of $`P`$ indicate that $`\eta 1`$, so $`B_{\mathrm{rnd}}`$ dominates the emissivity. To the extent that the random field is close to equipartition ($`ϵ_{B,\mathrm{rnd}}1`$), $`\eta =ϵ_{B,\mathrm{ord}}/ϵ_{B,\mathrm{rnd}}ϵ_{B,\mathrm{ord}}`$. If the shock is radiative during its early evolution, then cooling-induced compression increases the emissivity-weighted $`ϵ_{B,\mathrm{ord}}`$ over its immediate postshock (adiabatic) value by a factor888This assumes that the fraction $`ϵ_e`$ of the internal energy just behind the shock transition that resides in relativistic electrons and $`e^\pm `$ pairs is radiated away (Granot & Königl 2001). $`(1ϵ_e)^1`$. The transition from fast to slow cooling, which occurs at $`t=t_0`$, could therefore reduce $`\eta `$ and may contribute to the early decline of $`P`$ observed in some sources. During the subsequent, slow-cooling phase, $`ϵ_{B,\mathrm{ord}}`$ is essentially equal to the magnetization parameter of the ambient medium, $`\sigma =B_{\mathrm{ext}}^2/4\pi \rho _{\mathrm{ext}}c^2`$, so the evolution of $`P`$ during that phase may reflect the radial behavior of this parameter: $`\sigma `$ is expected to be roughly constant for an ISM or a stellar wind but to increase with $`r`$ inside a PWB (Königl & Granot 2002). If the orientation of the ambient field also changed with radius then this would lead to a gradual variation in $`\theta _p`$. If one approximates $`B_{\mathrm{ext}}r^{a/2}`$ and $`\rho _{\mathrm{ext}}r^k`$, then $`ϵ_{B,\mathrm{ord}}t_{\mathrm{obs}}^{\frac{a+k}{4k}}`$. We parameterize the above effects by $`ϵ_{B,\mathrm{ord}}=ϵ_{B,0}(t_{\mathrm{obs}}/t_0)^{\frac{a+k}{4k}}F(t_{\mathrm{obs}}/t_0)`$, where $`F(x)=1+\frac{g1}{g+1}\frac{2}{\pi }\mathrm{arctan}(\xi \mathrm{ln}x)`$ describes the amplitude ($`g`$) and sharpness ($`\xi `$) of the change in $`ϵ_{B,\mathrm{ord}}`$ at $`t_{\mathrm{obs}}t_0`$. We assume that $`ϵ_{B,\mathrm{tot}}=ϵ_{B,\mathrm{ord}}+ϵ_{B,\mathrm{rnd}}=\mathrm{const}`$, so $`\eta =(\frac{ϵ_{B,\mathrm{tot}}}{ϵ_{B,\mathrm{ord}}}1)^1`$.

Figure 3 shows an example of the polarization lightcurves. The choice of parameters for $`\eta (t)`$ is motivated by GRB 020813, in which $`\theta _p`$ was again roughly constant with time whereas $`P`$ first decreased (from $`P2\%`$ after $`6\mathrm{hr}`$ to $`P0.6\%`$ after $`24\mathrm{hr}`$) and subsequently increased monotonically, reaching $`P3.7\%`$ after $`96\mathrm{hr}`$. A sharp break in the lightcurve was observed after $`14\mathrm{hr}`$ (Covino et al. 2003b). It is seen that roughly equal contributions to the polarization from $`B_{\mathrm{ord}}`$ and $`B_{\mathrm{rnd}}`$ can provide a qualitative fit to the evolution of $`P`$ and $`\theta _p`$ in this source for $`\delta 60^{}90^{}`$. More generally, the polarization lightcurves can show a diverse behavior that varies as a function of $`\delta `$ as well as of $`b`$ and $`\theta _{\mathrm{obs}}/\theta _0`$. So long as $`\eta P_{\mathrm{ord}}>P_{\mathrm{rnd}}`$, the changes in $`\theta _p`$ would be small whereas the variations in $`P`$ could be significant, as found observationally.

## 4. Discussion

The linear polarization in GRB afterglows may be largely due to an ordered magnetic field in the ambient medium, which gives rise to an ordered field component behind the afterglow shock that is coherent over the entire emission region. This can result in a polarization position angle (PA), $`\theta _p`$, that is roughly constant in time as well as in a variable degree of polarization, $`P`$, as found in all afterglow observations to date (except one; see footnote 11).

The magnetic field in the GRB ejecta is potentially much more ordered than in the shocked ambient medium behind the afterglow shock, reflecting the likely presence of a dynamically important, predominantly tansverse, large-scale field advected from the source (e.g., Spruit, Daigne, & Drenkhahn 2001; Vlahakis & Königl 2001). This could result in a large value of $`P`$ \[up to $`(0.900.93)P_{\mathrm{max}}60\%`$\] in the prompt $`\gamma `$-ray emission999After this paper was submitted, Coburn & Boggs (2003) reported a measurement of $`P=80\pm 20\%`$ in the $`\gamma `$-ray emission of GRB 021206, which is naturally (and most likely; Granot 2003) produced in this way. as well as in the ‘optical flash’ and ‘radio flare’, which are attributed to emission from the reverse shock. If the polarization from the reverse shock is indeed dominated by the ordered component, and if it is coherent over the whole ejecta, then $`\theta _p`$ is not expected to vary significantly during the optical flash or between the optical flash and the radio flare. However, if the ordered magnetic field is coherent only in patches of angular size $`\theta _B`$, then, so long as $`\gamma >1/\theta _B`$, we expect $`PP_{\mathrm{max}}`$, whereas after $`\gamma `$ drops below $`1/\theta _B`$ we expect $`P\gamma \theta _BP_{\mathrm{max}}`$ and variations in $`\theta _p`$ on timescales $`\mathrm{\Delta }t_{\mathrm{obs}}t_{\mathrm{obs}}`$ on account of the averaging over $`N(\gamma \theta _B)^2`$ mutually incoherent patches within the observed region of angle $`1/\gamma `$ about the line of sight. (This resembles the proposal by Gruzinov & Waxman 1999, except that here $`N`$ is envisioned to increase with time.) In the latter case $`P`$ might be smaller and $`\theta _p`$ would be different in the ‘radio flare’ (for which typically $`\gamma 10`$) than in the ‘optical flash’ (for which $`\gamma 100`$).

Variability in the afterglow lightcurve, as reported in GRBs 021004 and 030329, whether induced by a clumpy external medium or a patchy shell (Lazzati et al. 2002; Nakar, Piran, & Granot 2002), should give a different weight to emission from different parts of the afterglow image, thus breaking its symmetry and inducing polarization.101010If the density distribution is spherically symmetric, the symmetry would need to be broken by the outflow geometry — e.g., a jet observed off axis. Therefore, we expect a highly variable lightcurve to be accompanied by variability in both $`P`$ and $`\theta _p`$.111111After this paper was submitted, a change of $`45^{}`$ in $`\theta _p`$ was reported in GRB 021004 between $`9`$ and $`16`$hr (Rol et al. 2003). This cannot be explained by simple jet models (Sari 1999; Ghisellini & Lazzati 1999) but could naturally arise in conjunction with the variability in the lightcurve (which, in fact, peaked at about the same time).

Early polarization measurements, starting at $`t_{\mathrm{obs}}t_j`$, are crucial for distinguishing between our model and purely tangled jet field models, as the latter predict $`P(t_{\mathrm{obs}}t_j)P(t_{\mathrm{obs}}t_j)`$, whereas our model allows $`P(t_{\mathrm{obs}}t_j)P(t_{\mathrm{obs}}t_j)`$. In the latter models $`P`$ is expected to peak, or else vanish and reappear rotated by $`90^{}`$, around $`t_j`$. In contrast, in our model, if the polarization is dominated by an ordered magnetic field, then the variations in the polarization around $`t_j`$ would be much less pronounced, with $`\theta _p`$ exhibiting only a gradual variation and $`P`$ never crossing zero. Our model predicts a possible change in $`P`$ around the transition time from fast to slow cooling, $`t_0`$, where typically $`t_01`$hr ($`1`$day) for ISM-like (stellar wind-like) parameters (although it may vary considerably around these values).

We thank P. Goldreich, R. Sari, A. Panaitescu, and E. Rossi for useful discussions. This research was supported in part by funds for natural sciences at the Institute for Advanced Study (JG) and by NASA ATP grant NAG5-12635 (AK).


format: text

text:
# Modelling the nova rate in galaxies

## 1 Introduction

Classical novae are binary systems, where one component is a white dwarf (either carbon-oxygen or neon-oxygen-magnesium), accreting hydrogen rich material from a less evolved companion. When the pressure and temperature at the bottom of the accreted layer exceed critical values (e.g. Livio 1992) the nuclear reactions are ignited and novae exhibit a sudden and rapid increase of their brightness. The thermonuclear runaways are responsible for the production of interesting chemical species, such as 13C, 15N, 17O, 22Na, 26Al, 22Ne and possibly 7Li (see Jose’ & Hernanz 1998, Gehrz et al. 1998, Della Valle et al. 2002) and for ejection of the accreted envelope. Although the matter processed through nova explosions is a small fraction (a few $`10^3`$) of the total mass of the interstellar gas and dust in the Milky Way, the concentrations of these rare isotopes in the nova ejecta can be enhanced, relative to solar abundances, by factors 102-3, therefore making novae non-negligible contributors to the galactic nucleosynthesis (Romano et al. 1999, 2001; Romano & Matteucci 2003). Novae have been observed not only in our Galaxy but also in the Local Group (e.g. Hubble 1929, Arp 1956 and Rosino 1964) and beyond, in Virgo (Pritchet & van den Bergh 1987) in M81 (Shara, Sandage & Zurek 1999) and in Fornax (Della Valle & Gilmozzi 2002). From these studies some estimates of nova rates have been derived and they are shown in Table 1, where the Hubble type of the galaxies is indicated in Col. 1, the nova rate in yr-1 in Col. 2, the K-band luminosity of the object in Col. 3 and in Col. 4 the luminosity specific nova rate in units of $`\mathrm{yr}^110^{10}\mathrm{L}_{,\mathrm{K}}`$ ($`LSNR_K`$). These rates are taken from Della Valle (2002) and Shafter et al. (2000).

In this paper we use the chemical evolution models elaborated by Chiappini et al. (2001), Matteucci et al. (1998) and Calura et al. (2003) for the Milky Way, Ellipticals and late Spirals, respectively, to compute the nova rate in ‘disk’ dominated and ‘bulge’ dominated systems with the intent of checking whether the different observed nova rates can be reproduced under the assumption that the main difference between galactic morphological types should be ascribed to the different histories of star formation. Moreover, our study is aimed at deriving important constraints on some key parameters of nova theory, which cannot be directly inferred from observations, such as the recurrence time between the outbursts, the fraction of binary systems which give rise to nova systems, the delay time with which the first novae appear relative to the beginning of star formation. The use of such models is appropriate since they can reproduce the basic parameters characterizing the respective parent galaxies such as the rate of type Ia SNe, B and K luminosities of the parent galaxies, the \[$`\alpha `$/Fe\] ratios and global metallicities in the stellar populations and gas. Finally, we note that a fraction of type Ia SNe can originate from WDs in CV-like systems, after reaching the Chandrasekhar mass limit (e.g. Livio 2000). Whether these systems can give rise to novae or supernovae depends crucially on the mass of the underlying white dwarf and mass accretion rate (see Nomoto et al., 1984). Therefore, when trying to predict nova rates in galaxies it is obvious to compute and try to reproduce the observed SNIa rates. The observed type Ia rates in galaxies of different morphological type are shown in Table 2.

While the rates of novae are expressed in $`\mathrm{yr}^1`$ or in $`LSNR_K`$, the rates of SNIa are expressed in $`SNu_B`$ (SN $`(100\mathrm{y}\mathrm{r})^110^{10}\mathrm{L}_\mathrm{B}_{}`$) and $`SNu_K`$, (SN $`(100\mathrm{y}\mathrm{r})^110^{10}\mathrm{L}_\mathrm{K}_{}`$). The paper is organized as follows: in Sect. 2 we describe how to compute the nova and SN rates in galaxies and compare the results with observations and finally, in Sect. 3, we draw some conclusions.

## 2 Modelling the nova and supernova rates

### 2.1 The Milky Way

The model for the Galaxy assumes two main infall episodes for the formation of the halo-thick disk, and the thin-disk, respectively. The timescale for the formation of the thin-disk is much longer than that of the halo, implying that the infalling gas forming the thin-disk comes not only from the halo but mainly from the intergalactic medium. The timescale for the formation of the thin disk is assumed to be a function of the galactocentric distance, leading to an inside-out picture for the Galaxy disk build-up, according to the original suggestion of Matteucci & François (1989). The two-infall model differs from other models in the literature mainly in two aspects: it considers an almost independent evolution between the halo and thin-disk components (see also Pagel & Tautvaisiene 1995), and it assumes a threshold in the star formation process, (see Kennicutt 1989, 1998; van der Hulst et al. 1993; Martin & Kennicutt 2001). The model well reproduces the majority of observational constraints about the abundances of heavy elements both locally and in the whole disk.

The star formation rate (SFR) adopted here has the same formulation as in Chiappini et al. (1997):

$`\psi (R,t)=\nu \left(t\right)\left({\displaystyle \frac{\mathrm{\Sigma }(R,t)}{\mathrm{\Sigma }(R_{},t)}}\right)^{2\left(k1\right)}\left({\displaystyle \frac{\mathrm{\Sigma }(R,t_{Gal})}{\mathrm{\Sigma }(R,t)}}\right)^{k1}G_{gas}^k(R,t)`$ (1)
where $`\nu (t)`$ is the efficiency of the star formation process, $`\mathrm{\Sigma }(R,t)`$ is the total surface mass density at a given radius $`R`$ and given time $`t`$, $`\mathrm{\Sigma }(R_{},t)`$ is the total surface mass density at the solar position, $`G_{gas}(R,t)`$ is the surface gas density normalized to the present time total surface mass density in the disk $`\mathrm{\Sigma }_D(R,t_{Gal})`$, tGal = 13 Gyr is the age of the Galaxy, R = 8 kpc is the solar galactocentric distance (see Reid 1993).

The gas surface density exponent, $`k`$, is set equal to 1.5, in order to ensure a good fit to the observational constraints in the solar vicinity. This value is also in agreement with the observational results of Kennicutt (1998), and with N-body simulation results by Gerritsen & Icke (1997). The star formation efficiency is set to $`\nu `$ = 2 Gyr-1 for the Galactic halo, whereas it is $`\nu `$ = 1 Gyr-1 for the disk; this is to ensure the best fit to the observational features in the solar vicinity. The star formation rate becomes zero when the gas surface density drops below a certain critical threshold. We adopt such a threshold density to be $`(\sigma _{gas})_{th}`$ $``$ 4 $`M_{}`$ pc-2 for the Galactic halo, and $`(\sigma _{gas})_{th}`$ $``$ 7 $`M_{}`$ pc-2 for the disk (see Chiappini et al. 2001). The assumption of such threshold densities naturally produces the existence of a hiatus in the SFR between the halo-thick disk phase and the thin-disk phase. This discontinuity in the SFR is observed in the \[Fe/O\] vs. \[O/H\] (Gratton et al. 2000) and in the \[Fe/Mg\] vs. \[Mg/H\] (Fuhrmann 1998) plots. The IMF is that of Scalo (1986), assumed to be constant in time and space. We compute the nova systems formation rate ($`R_{novae}`$(t)) at the time $`t`$ as a fraction of the formation rate of white dwarfs at a previous time $`t\mathrm{\Delta }t`$ (see D’Antona & Matteucci, 1991):

$`R_{novae}(t)=\alpha {\displaystyle _{0.8}^{8.0}}\psi (t\tau _m\mathrm{\Delta }t)\varphi (m)𝑑m.`$ (2)
Here $`\mathrm{\Delta }t`$ is a delay time whose value has to be fixed to guarantee the cooling of the white dwarf (WD) to a level that ensures a strong enough nova outburst. We assume a distribution of $`\mathrm{\Delta }t`$ varying from 1 to 5 Gyr. It is worth noting that assuming a constant average delay of $`\mathrm{\Delta }t`$=2 Gyr, as assumed in previous papers (e.g. Romano et al. 1999; 2001), produces negligible differences in the final results.

It is assumed that all stars with initial masses in the range 0.8 - 8 M end their lives as WDs. $`\psi `$(t) is the SFR as defined in Eq. (1), $`\tau _m`$ is the lifetime of the star of mass m and $`\varphi (m)`$ is the IMF.

The rate of nova eruptions is related to the WD formation rate by:

$`R_{outbursts}(t)=\alpha R_{WDs}(t)n,`$ (3)
where $`\alpha R_{WDs}=R_{novae}(t)`$ is the formation rate of WDs in binary systems, as defined in Eq. (2), which will give rise to nova eruptions and $`n=10^4`$ is the average number of nova outbursts for a typical nova system (Ford, 1978; Bath & Shaviv 1978; Shara et al. 1986). This number, which is just an average value, is inversely proportional to the recurrence time. The parameter $`\alpha `$, which here is a constant in space and time, can be fixed by reproducing the rate of nova outbursts in our Galaxy at the present time. Unfortunately, estimates of this quantity in the current literature show a large spread. Predictions based on scalings from extragalactic nova surveys suggest low values, whereas estimates based on extrapolations of Galactic nova observations give the highest rates. We consider for the present time observed rate of nova outbursts in the Galaxy $`R_{outbursts}(t_{Gal})25\mathrm{y}\mathrm{r}^1`$, as in Romano et al. (1999), for the following reasons. Given the fact that observations of novae in nearby galaxies would avoid, or at least minimize, most of the problems encountered by Galactic observations, such as interstellar extinction in the Galactic disk. Shafter (1997) shows that the nova rate based on Galactic observations can be made consistent with the rate predicted from the extragalactic data (Della Valle et al. 1994) particularly if the Galaxy has a strong bar oriented in the direction of the Sun. In this latter, most favourable case, the suggested value is near $`20\mathrm{y}\mathrm{r}^1`$, otherwise, if the bar is weak or misaligned, the global rate can be $`30\mathrm{y}\mathrm{r}^1`$.

Our first choice is to fix $`\alpha `$ = 0.01 (coupled with $`n=10^4`$) so to reproduce the present observed nova rate of 25$`\mathrm{yr}^1`$ and to use the same $`\alpha `$ in the models for elliptical and irregular galaxies. This value of $`\alpha `$ is the same as in previous papers (e.g. Romano et al. 2001; Matteucci et al. 1999).

Finally, the SNeIa rate has been computed following Greggio & Renzini (1983) (see also Matteucci & Greggio 1986) and is expressed as:

$`R_{SNeIa}=A{\displaystyle _{M_{B_m}}^{M_{B_M}}}\varphi (M_B){\displaystyle _{\mu _m}^{0.5}}f(\mu )\psi (t\tau _{M_2})𝑑\mu 𝑑M_B,`$ (4)
where $`M_2`$ is the mass of the secondary, $`M_B`$ is the total mass of the binary system, $`\mu =M_2/M_B`$, $`\mu _m=max\{M_2(t)/M_B,(M_B0.5M_{B_M})/M_B\}`$, $`M_{B_m}`$ = 3 $`\mathrm{M}_{}`$, $`M_{B_M}`$ = 16 $`\mathrm{M}_{}`$. $`\varphi (M_B)`$ is the initial mass function for the total mass of the binary system, $`f(\mu )`$ is the distribution function for the mass fraction of the secondary, $`f(\mu )=2^{1+\gamma }(1+\gamma )\mu ^\gamma `$, with $`\gamma `$ = 2; $`A`$ = 0.05 is the fraction of systems with total mass in the appropriate range, which give rise to SNIa events. This quantity is fixed by reproducing the observed SNe Ia rate at the present epoch (see also Madau et al. 1998), and it provides excellent fits to the abundances and abundance ratios in the Galaxy.

The results for the Galaxy are shown in Table 3 and Fig. 1. In Table 3 we report the predicted nova and SN Ia rate for the Milky Way, adopting an age of 13 Gyr, according to a cosmology with $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$ and $`z_f=10`$. We computed the K- and B- band luminosity of the Milky Way by means of the spectro-photometric model of Jimenez et al. (1999): we find, for the present time, luminosities in very good agreement with the estimated ones. The estimated Galactic blue luminosity is $`\mathrm{L}_{\mathrm{B}_{\mathrm{Gal}}}210^{10}\mathrm{L}_{}`$ (see van den Bergh 1988) and the K-luminosity is $`\mathrm{L}_{\mathrm{K}_{\mathrm{GAL}}}1210^{10}\mathrm{L}_{}`$ (see Shafter el al. 2000). In Figure 1 we show the behaviour of the star formation rate in the framework of the two-infall model, where the gap between the formation of halo-thick disk and thin-disk is evident, as well as the behaviour of the nova and SN Ia rates. These latter show only a small discontinuity in correspondance with the star formation gap and the reason is that they both depend upon the past star formation. The nova rate increases continuously up to the present time due to the long delay assumed for these systems (the lifetime of the WD plus the $`\mathrm{\Delta }t`$), whereas the SNIa rate reaches a maximum and then it stays flat afterwards.

### 2.2 Ellipticals

In the models for ellipticals, the galaxies can be considered initially as an homogeneous sphere of gas with luminous masses in the range $`10^910^{12}\mathrm{M}_{}`$. A single zone interstellar medium with instantaneous mixing of gas is assumed throughout.

The adopted model parameters are shown in Tables 4 and 5, where we show the initial luminous masses for the model galaxies, the assumed star formation efficiency parameter, the effective radius, the predicted time for the occurrence of a galactic wind $`t_{GW}`$ and the predicted K- and B-luminosities. The parameter $`A`$ is varied according to the assumed IMF in order to reproduce the present time observed type Ia SN rate (Matteucci & Gibson, 1995). On the other hand, we kept the parameter $`\alpha `$ to be the same as in the Milky Way. The model is described in Matteucci et al. (1998) where we address the reader for more details.

The star formation rate is given by:

$`\psi (t)=\nu \rho _{gas}(t)/\rho (0),`$ (5)
i.e. normalized to the initial total volume density. $`\psi (t)`$ is assumed to drop to 0 at the onset of the galactic wind. The quantity $`\nu `$ is expressed in units of $`\mathrm{Gyr}^1`$ and represents the efficiency of star formation, namely the inverse of the time scale of star formation. The star formation is assumed to stop after the development of a galactic wind occurring before than 1 Gyr, from the beginning of star formation, for all the galaxies listed above. Therefore, the star formation rate in these galaxies can be considered as a strong burst which does not last more than 1 Gyr and is shorter in more massive systems. This is obtained by assuming that the star formation efficiency increases with galactic mass thus producing an “ inverse wind” effect, as described in Matteucci (1994), where the galactic wind occurs before in more massive than in less massive ellipticals. As a consequence of this, the star formation period is longer in smaller systems thus allowing the SNe Ia to substantially pollute the ISM. This effect can explain the observed increase of the \[$`\alpha `$/Fe\] ratio with galactic mass (Worthey et al. 1992; Matteucci 1994), which is not obtained in a classic wind scenario, where the more massive objects form stars for a longer time (Larson, 1974). The galactic wind develops as a consequence of the energy transfer from SNe into the ISM. In fact, when the thermal energy of the gas becomes larger than the binding energy of the gas, the wind starts (Arimoto & Yoshii, 1987; Matteucci & Tornambe’ 1987; Matteucci 1992, 1994; Pipino et al. 2002). In order to compute the binding energy of the gas some assumptions have to be made about the galactic potential well. In particular, it is assumed that all ellipticals possess heavy but diffuse dark matter halos; a ratio between the half-light radius and the radius of the dark matter core $`R_{luminous}/R_{dark}`$=0.10 and a ratio dark to luminous mass of 10 are assumed.

The B and K luminosities for the ellipticals of different masses are computed by means of the spectro-photometric model of Jimenez et al. (1999) and are used to compute the SN Ia rate in SNu and the nova rate per unit of $`L_K`$, respectively.

The main differences between the model for the Galaxy, and the model for an elliptical galaxy concern the different SFR, which is much stronger in the earliest stages, and then is set to zero after the galactic wind in the case of ellipticals. It is worth noting that galactic winds in the case of elliptical galaxies seem necessary to explain their lack of gas and the chemical enrichment of the ICM, whereas the Galactic model does not take into account the occurrence of a strong wind, mainly because of the strong gravitational potential well associated with our Galaxy, and the presence of the gas in the Galactic disk.

We adopt two different initial mass functions (IMFs): the Salpeter ($`x_{IMF}=1.35`$) one and the Arimoto & Yoshii (1987) ($`x_{IMF}=0.95`$) one. In fact, successful models of chemical evolution of ellipticals have shown that the Scalo IMF is not suitable for these galaxies (e.g. Matteucci & Gibson, 1995).
Table 6 shows the results of the model I (Salpeter IMF), whereas Table 7 shows the results of the model II (x=0.95 IMF). These models have been computed by assuming the same $`\alpha `$ and the same $`n`$ as for the Milky Way.

In particular, in the first column we report the initial galactic luminous mass, in the second column the nova rate in units of $`\mathrm{yr}^1`$, in Col. 3 the luminosity specific nova rate , in Col. 4 the SN Ia rate in units of $`(100\mathrm{y}\mathrm{r})^1`$ and in Col. 5 and 6 the SN Ia rate in units of $`SNu_B`$ and $`SNu_K`$, respectively. The predicted nova rates are quite large for massive ellipticals ranging between 300 and 1000 nova $`\mathrm{yr}^1`$ which is about a factor 3-10 larger than it is derived from observations. If these values were realistic, the discrepancy could be partially due to an observational bias affecting the ground based nova surveys due to poor spatial resolution and bright limiting magnitude. An indication in this direction comes from the nova rate of M87 recently provided by Shara & Zurek (2002) on the basis of HST archive images, which is a factor 2–3 larger than previous ground-based estimates (see Table 1). Another possibility is that in early type galaxies the recurrence time between two consecutive nova explosions is considerably longer than in late spirals as a consequence of the different stellar population from which novae originate. This suggestion is supported both by observational and theoretical grounds. From one side Duerbeck (1990) and Della Valle et al. (1992, 1994) have demonstrated, on the basis of galactic and extragalactic nova observations, the existence of two populations of novae: fast and bright novae belonging to ‘disk’ stellar population and slow and faint novae which originate from a ‘bulge’ stellar population. In particular, the latter authors (see also Della Valle & Livio 1998) suggested that bulge novae could originate from relatively light WDs, likely in the range of masses of $`M_{WD}0.9M_{}`$ while novae in the disk arise from massive WDs ($`M_{WD}1M_{}`$). From the other side, Truran (1990, see also Ritter et al. 1991) has found that the mass of the WD and the recurrence time between the outbursts are inversely proportional. In order to match the ‘theoretical’ with ‘empirical’ rates, one needs to lower either $`\alpha `$ or $`n`$ or both. For example for M87 the match between predicted and observed rates can be achieved by increasing the recurrence time between two consecutive outbursts, to $`T_r3110^5`$ yr. This is about 10–3 times larger than assumed for novae in the ‘disk’ of the Galaxy. In this way, both the predictions of models I and II, for a luminous mass in the range 4-6 $`10^{11}M_{}`$ are in very good agreement with the preliminary new estimate from Shara & Zurek (2002) of the nova rate for M87. In fact, M87 can be modeled both as an elliptical galaxy with luminous mass $`410^{11}M_{}`$ and $`R_{eff}6`$ kpc or with luminous mass $`610^{11}M_{}`$ and $`R_{eff}`$ 7 kpc (Cohen & Ryzhov 1997). The predictions of these two models do not differ much and both of them can well represent M87. Another possibility is to adopt a fraction of binary systems giving rise to novae smaller than 0.01. Unfortunately, if from one side we can justify this value for the Milky Way because this $`\alpha `$ gives both the observed present time Galactic rate and the correct yields of Lithium and CNO isotopes (see Romano & Matteucci, 2003), on the other side, we completely lack this information for ellipticals.

Finally, a comparison between tables 6 and 7 with Table 2 shows a very good agreement between the observed and predicted type Ia SN rates both in units of $`SNu_B`$ and $`SNu_K`$.

In figures 2 and 3 we show the predicted star formation, nova and SNIa rates as functions of the cosmic time for ellipticals of different luminous masses and for the two IMFs.

Figure 4 shows the nova rate in the ellipticals as a function of the mass of the galaxy. From this figure one can infer the expected nova rate as a function of luminous galactic mass.

### 2.3 Magellanic Irregulars

As an example of Magellanic irregular galaxies and late-type spiral we computed a model to fit the LMC. In particular, we assumed an history of star formation with two main bursts occurring during the first 3 Gyr of the galactic life and during the last 2 Gyr, respectively, coupled with a lower star formation rate in between the two, as shown in Figure 5. This star formation regime for the LMC has been suggested by Calura et al. (2003) and it can reproduce the observed \[O/Fe\] versus \[Fe/H\]. The star formation rate is like in Eq. (5) with an efficiency of star formation $`\nu =0.1\mathrm{Gyr}^1`$ in each starburst and of $`\nu =0.01\mathrm{Gyr}^1`$ in the interburst period. The IMF is the Salpeter one and the parameters $`\alpha =0.01`$ and $`A=0.18`$. The history of star formation has been chosen to reproduce \[O/Fe\] vs. \[Fe/H\] in the LMC (see Figure 5 and Calura et al. 2003). We considered a luminous mass of $`10^{10}M_{}`$ which is suitable for the LMC (e.g. Russell & Dopita, 1992). We adopted the same values of $`\alpha `$ and $`n`$ as for the Milky Way and found that the predicted present time nova rate is in excellent agreement with the observed one and is $`1.7`$ nova $`\mathrm{yr}^1`$. The predicted B and K luminosities are $`L_B=210^9L_{}`$ and $`L_K=0.2710^{10}L_{}`$, respectively. As a consequence, we predict a nova rate in units of $`LNSR_K`$ of $`6`$, in agreement with the data of Table 1. The predicted type Ia SN rate is $`0.31`$SNu, also in very good agreement with the data of Table 2.

## 3 Discussion and Conclusions

In this paper we have computed the Galactic and extragalactic nova rates, in particular the nova rates expected in elliptical galaxies and in Magellanic irregulars, typified by the LMC. To do that we have adopted detailed models of chemical and photometric evolution of the Milky Way, ellipticals and the LMC. These models can reproduce the majority of the observational constraints, including abundances, abundance ratios, luminosities and colors. We have computed also the type Ia SN rates since the novae and type Ia SNe are originating by CV-like binary systems. The nova rate is calculated in a simple way but still depends on several unknown parameters such as the fraction of binary systems which are nova progenitors (parameter $`\alpha `$) and the total number of outbursts occurring during a nova lifetime (parameter $`n`$). For the sake of simplicity, we have assumed these two parameters to be the same in all galaxies whereas the star formation histories are different and tuned to reproduce the available observational constraints. For the Milky Way this particular choice of parameters, besides reproducing the observed nova rate, allows to predict $`{}_{}{}^{7}Li`$, $`{}_{}{}^{15}N`$, $`{}_{}{}^{13}C`$ and $`{}_{}{}^{17}O`$, in good agreement with the observations (Romano & Matteucci, 2003). On the other hand, for external galaxies we do not have other constraints related to novae except the nova rates. Therefore, this approach does not allow us to put strong constraints upon $`\alpha `$ and $`n`$ singularly but rather on the product $`\alpha n`$. Concerning the delay with which novae appear since the beginning of star formation, we find that an average delay of 2 Gyr can be assumed for all the studied galaxies.

Our main conclusions can be summarized as follows:

* We predict that the nova rates in ellipticals reach a maximum after $`10^9`$ years since the beginning of star formation and then they decrease after that. The same behaviour is shared by type Ia SNe although they reach a maximum after only a few $`10^8`$ years (see also Matteucci & Recchi, 2001). The delayed maximum in the nova rate relative to the SNIa rate is due to the longer time delay assumed for novae to appear relative to SNeIa (see Canal et al. 1996).
* Our results are suggesting that a model for ellipticals which assumes a strong burst of star formation occurring at early cosmic times can reproduce the majority of the observational constraints in ellipticals. The same model predicts that the nova rates are increasing with the galactic mass of ellipticals and for a galaxy like M87 we predict a present time nova rate of 300-1000 nova $`\mathrm{yr}^1`$, when the values of $`\alpha `$ and $`n`$ are the same as for the Galaxy. A good agreement with the observations (100-300 nova $`\mathrm{yr}^1`$) is obtained if we assume an average recurrence time between two consecutive nova outbursts of the order of $`1310^5`$ yr, which is a factor of 10 larger than the typical recurrence time for novae originating in the disk of the Milky Way. This fact is a consequence of the existence of two populations of novae, one associated with disk/spiral arms stellar population and the other with bulge/thick disk stellar population (Della Valle et al., 1994; Yungelson et al. 1997).
* The nova rate in ellipticals is rather independent from the adopted IMF, but it depends on the assumption that the efficiency of star formation is an increasing function of the luminous galactic mass (inverse wind model of Matteucci 1994). We do not predict any trend for the luminosity specific nova rate as a function of the galactic mass for these galaxies.
* We predict also a type Ia SN rate for ellipticals in very good agreement with the observed estimates and a slight decreasing trend for the SNIa rate with galactic blue and red luminosity.
* Finally, we confirm the existence of a trend for the luminosity specific nova rate, as derived from the observations (Della Valle et al. 1994), passing from a Magellanic irregular (spiral) like the LMC, $`5`$, to a Sb spiral galaxy like the Milky Way $`2.1`$. For a giant elliptical like M87 we obtain a large range of uncertainty between 2.5 and 7.5 novae/yr per unit of 1010K luminosity. The lower limit of the rate is very similar to that obtained from ground based nova survey whereas the upper limit is consistent with nova rates obtained from HST observations (Shara & Zurek 2002). If true, this could indicate that the studies of novae in extragalactic systems from ground-based telescopes, can be largely affected by observational bias. As an alternative, we note that Livio et al. (2002) have suggested that M87 could be a fast producer of novae, because of the presence of the jet which can enhance the accretion onto white dwarfs and therefore increase the nova rate.

###### Acknowledgements.

We thank Donatella Romano for many useful suggestions and the referee H.W. Duerbeck for his careful reading of this manuscript.


text:
# On exact polytropic equilibria of self-gravitating gaseous and radiative systems: their application to molecular cloud condensation

## 1 Introduction

Fundamental research on the structure and evolution of stars, molecular clouds, and galaxies is still very important (Longair, 1994). As a rule, a family of the Lane-Emden (LE) equations classified by the polytrope indices $`n`$ is known to be the most powerful tool to investigate the self-gravitating polytropic equilibria (Binney & Tremaine, 1987). So far, many numerical and semi-analytical works have been devoted to seeking solutions (Liu, 1996; Roxburgh & Stockman, 1999; Goenner & Havas, 2000; Hunter, 2001; Medvedev & Rybicki, 2001), which include the non-singular function with $`n\mathrm{}`$ for an isothermal polytropic sphere (e.g., Natarajan & Lynden-Bell, 1997). Removing infinite divergence of the enclosed mass, some modified models are used to depict the globular clusters, intergalactic media, clusters of galaxies, and so forth (King, 1962; Jones & Forman, 1984).

Prior to modern numerical methods (VandenBerg & Bell, 1985; VandenBerg, 1985), the pioneering stellar model proposed by Eddington, which is based on the solution of the LE equation with the polytrope index of $`n=3`$, could be useful to envisage stellar interiors, although the complexities such as transport processes and non-ideal equation of state (EOS) have been omitted. Regarding the study of the interstellar medium (ISM), it seems that some classes of the LE equations shed light on dark molecular clouds and their condensations (Mizuno et al., 1994). However, in the general case that the pressure which supports the celestial objects is provided largely by the thermal motions of particles, we need to pay attention to the extensive use of the LE equation with a specified index. For heuristic ways, we explain this crucial point below.

Let us suppose that the object sustains its own self-gravity due to the total pressure of $`P=p_\mathrm{g}+p_\mathrm{r}`$, where the partial pressures of gas and radiation field are described by the EOS of $`p_\mathrm{g}=\beta P=[k_\mathrm{B}/(\mu m_\mathrm{p})](T/V)`$ and $`p_\mathrm{r}=(1\beta )P=\frac{1}{3}aT^4`$, respectively. Here, $`m_p`$ is the proton mass, $`V=\rho ^1`$ is the specific volume, $`a=8\pi ^5k_\mathrm{B}^4/(15c^3h^3)`$ is the Stefan-Boltzmann constant, and $`\mu 1`$ denotes the mean molecular weight. 111For plasmas, $`1/2\mu 2/(1+3X+0.5Y)2`$, where the mass fractions contained in the unit mass of the medium are $`X\mathrm{g}`$ hydrogen, $`Y\mathrm{g}`$ helium, and $`Z=1XY\mathrm{g}`$ heavy elements. For $`0<\beta <1`$, the relation between temperature and density is given by $`T=\{[k_\mathrm{B}/(\mu m_\mathrm{p})](3/a)[(1\beta )/\beta ]\}^{1/3}V^{1/3}`$. Hence, we obtain the polytropic relation as follows:

$$P(r)=\left[\left(\frac{k_\mathrm{B}}{\mu m_\mathrm{p}}\right)^4\frac{3}{a}\frac{1\beta }{\beta ^4}\right]^{1/3}V^{4/3}(r).$$
(1)
If we assume $`\beta =\mathrm{const}`$ (Chandrasekhar, 1967), equation (1) has the form of $`P=K\rho ^{4/3}`$, where $`K=\{[k_\mathrm{B}/(\mu m_\mathrm{p})]^4(3/a)[(1\beta )/\beta ^4]\}^{1/3}`$ is the adiabatic constant. In this aspect, such a radiative polytrope leads to the LE equation of index $`n=3`$ for gravitational equilibria. Hereafter, we refer to this scheme, namely, the Eddington model, as ’LE$`3`$’, without further notice. Apparently, for gaseous objects of $`\beta 1`$, equation (1) does not match with the relation of $`P(r)\rho ^{5/3}(r)`$ involving $`n=3/2`$ for a perfect monatomic gas.

The possible scenario that recovers the physical consistency is to allow the spatial variation of polytrope index. In equation (1) such changes can be achieved by invoking the density dependence of $`\beta =\beta [V(r)]`$, and a trivial solution of $`\mathrm{d}\mu /\mu =0`$ (not shown). Of course, these are not mathematical artefacts. In the astronomical context, the present scheme is applicable to study, especially, the high-density adiabatic cores deeply immersed in the isothermal molecular cloud, which are so-called ’first cores’ (Masunaga, Miyama & Inutsuka, 1998). As yet, there has been no distinct detection of the first cores. We hope for future-planed observations of molecular cloud cores in the range of millimetre and submillimetre wavelengths with high angular resolution (Hasegawa, 2001).

In this paper, we present a self-consistent solution of a polytropic sphere for self-gravitating gaseous and radiative systems, going beyond the LE framework. In the following, we derive a generic equation for hydrostatic equilibria. In the present scheme, if we know the whole mass and the mean molecular weight of the system, then the pressure ratio of gas/total pressure at the centre, which contributes to the boundary condition, is identified. It is shown that the polytrope varies radially, modifying the density and pressure profiles derived from the LE$`3`$ model. In the limit that the gaseous pressure is asymptotically close to the total one, the solution exactly agrees with the LE solution of index $`n=3/2`$.

## 2 The nonlinear Poisson equation
including a self-consistent variable
polytrope

We begin by defining the adiabatic exponent as $`\mathrm{\Gamma }_1=\mathrm{d}(\mathrm{ln}P)/\mathrm{d}(\mathrm{ln}\rho )`$. Introducing the total derivative of $`\mathrm{d}P=\mathrm{d}(p_\mathrm{g}+p_\mathrm{r})=p_\mathrm{g}(\mathrm{d}V/V)+(p_\mathrm{g}+4p_\mathrm{r})(\mathrm{d}T/T)`$, we obtain

$$(p_\mathrm{g}+4p_\mathrm{r})\frac{\mathrm{d}T}{T}[p_\mathrm{g}(p_\mathrm{g}+p_\mathrm{r})\mathrm{\Gamma }_1]\frac{\mathrm{d}V}{V}=0.$$
(2)
For a quasi-static and adiabatic change, the thermodynamical principle requires

$$\frac{\mathrm{d}Q}{V}=\left(\frac{1}{\gamma 1}p_\mathrm{g}+12p_\mathrm{r}\right)\frac{\mathrm{d}T}{T}+(p_\mathrm{g}+4p_\mathrm{r})\frac{\mathrm{d}V}{V}0,$$
(3)
where $`\gamma =c_\mathrm{p}/c_\mathrm{v}`$ is the specific-heat ratio. Combining equation (2) with equation (3), we find the relation between $`\mathrm{\Gamma }_1`$ and $`\beta `$ in the form of

$$\mathrm{\Gamma }_1[\beta (V),\gamma ]=\beta (V)+\frac{\left(\gamma 1\right)\left[43\beta (V)\right]^2}{\beta (V)+12\left(\gamma 1\right)\left[1\beta (V)\right]}.$$
(4)
It should be noticed that equation (4) holds the asymptotic properties of $`\mathrm{\Gamma }_14/3`$ and $`\gamma `$ for $`\beta 0`$ and $`1`$, respectively.

On the other hand, taking the derivative of equation (1), we obtain $`\mathrm{d}P/\mathrm{d}V=(P/3)\{(43\beta )/[\beta (1\beta )]\}(\mathrm{d}\beta /\mathrm{d}V)(4/3)(P/V)`$. This may be rewritten as

$$\frac{\mathrm{d}P}{P}+\frac{4}{3}\frac{\mathrm{d}V}{V}+\frac{43\beta }{3\left(1\beta \right)}\frac{\mathrm{d}\beta }{\beta }=0.$$
(5)
According to the relation of $`\mathrm{d}P/P+\mathrm{\Gamma }_1\mathrm{d}V/V=0`$, therefore, we obtain the new relation of

$$V\frac{\mathrm{d}\beta }{\mathrm{d}V}=\frac{3\beta (1\beta )}{43\beta }\left[\mathrm{\Gamma }_1(\beta ,\gamma )\frac{4}{3}\right],$$
(6)
which effectively changes the polytrope index of equation (1).

Now we consider a hydrostatic equilibrium of a self-gravitating system, described by the spherically symmetric Poisson equation, namely, $`\mathrm{\Phi }=r^2(\mathrm{d}/\mathrm{d}r)(r^2\mathrm{\Phi })=4\pi G\rho `$, where $`\mathrm{\Phi }(r)=\rho ^1(\mathrm{d}P/\mathrm{d}r)`$ and $`G`$ is the gravitational constant. By using equations (1) and (6), the non-linear equation can be cast to

$$\frac{1}{4\pi G}\left[\left(\frac{k_\mathrm{B}}{\mu m_\mathrm{p}}\right)^4\frac{3}{a}\right]^{1/3}$$

$$\times \frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2\left[\frac{1\beta (\rho )}{\beta ^4(\rho )}\right]^{1/3}\frac{\mathrm{\Gamma }_1[\beta (\rho ),\gamma ]}{\rho ^{2/3}(r)}\frac{\mathrm{d}\rho (r)}{\mathrm{d}r}\right)=\rho (r).$$
(7)
In order to normalize equation (7), we introduce $`\stackrel{~}{\rho }(\xi )=\rho (\xi )/\rho _\mathrm{c}`$, where $`\rho _\mathrm{c}`$ is the central mass density, and we properly choose the dimensionless frame of $`\xi =r/\alpha `$. Here, the characteristic scalelength $`\alpha `$ is given by

$$\alpha (\rho _\mathrm{c},\mu )=\left[\frac{1}{(4\pi G)^3\rho _\mathrm{c}^2}\left(\frac{k_\mathrm{B}}{\mu m_\mathrm{p}}\right)^4\frac{3}{a}\right]^{1/6}.$$
(8)
Then, we newly find a set of ordinary differential equations in the following form:

$$\frac{1}{\xi ^2}\frac{\mathrm{d}\stackrel{~}{m}(\xi )}{\mathrm{d}\xi }=\stackrel{~}{\rho }(\xi ),$$

$$\stackrel{~}{m}(\xi )=\frac{\xi ^2}{\stackrel{~}{\rho }^{2/3}(\xi )}\frac{\mathrm{d}\stackrel{~}{\rho }(\xi )}{\mathrm{d}\xi }\left[\frac{1\beta (\xi )}{\beta ^4(\xi )}\right]^{1/3}\mathrm{\Gamma }_1[\beta (\xi ),\gamma ],$$

$$\frac{1}{\stackrel{~}{\rho }(\xi )}\frac{\mathrm{d}\stackrel{~}{\rho }(\xi )}{\mathrm{d}\xi }=\frac{43\beta (\xi )}{3\beta (\xi )[1\beta (\xi )]\{\mathrm{\Gamma }_1[\beta (\xi ),\gamma ]\frac{4}{3}\}}\frac{\mathrm{d}\beta (\xi )}{\mathrm{d}\xi }.$$
(9)
Note that the enclosed mass defined by $`m(r)=4\pi _0^r\rho (r^{})r^2dr^{}`$ is transformed into $`m(\xi )=4\pi \rho _\mathrm{c}\alpha ^3(\rho _\mathrm{c},\mu )\stackrel{~}{m}(\xi )0`$. In the special case of $`\mathrm{d}\beta /\mathrm{d}\stackrel{~}{\rho }0`$, equation (9) reduces to the LE equation with an index of $`n(\gamma 1)^1`$ and $`3`$ for $`\beta 1`$ and $`0`$, respectively.

## 3 Numerical solutions of the generalized equation

To survey the inner structure of the polytropic sphere, we numerically integrate equation (9) for $`\gamma =5/3`$, and the boundary conditions of $`\stackrel{~}{\rho }(\xi )|_{\xi 0}=1`$, $`\stackrel{~}{m}(\xi )|_{\xi 0}=0`$, and $`\beta (\xi )|_{\xi 0}=\beta _\mathrm{c}`$ where $`\beta _\mathrm{c}=(p_\mathrm{g}/P)_\mathrm{c}`$ is a parameter given at the centre. It is noted that the solutions of $`\stackrel{~}{\rho }(\xi )`$ and $`\stackrel{~}{m}(\xi )`$ tend to monotonically decrease outwards, while $`\beta (\xi )`$ monotonically increases. At $`\xi =\xi _1`$, $`\stackrel{~}{\rho }(\xi )`$ and $`\mathrm{d}\stackrel{~}{m}(\xi )/\mathrm{d}\xi `$ vanish, indicating a well-defined radius and mass of the object, i.e.

$$R=\alpha (\rho _\mathrm{c},\mu )\xi _1(\beta _\mathrm{c}),$$
(10)
$$M=m[\xi _1(\beta _\mathrm{c})]=4\pi \rho _\mathrm{c}\alpha ^3(\rho _\mathrm{c},\mu )\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})].$$
(11)
Table 1 lists the typical values of $`\xi _1`$ and $`\stackrel{~}{m}(\xi _1)`$ as a function of $`\beta _\mathrm{c}`$.

### 3.1 Inner structure of gaseous and radiative polytropic sphere

Substituting equation (8) into equation (11), we obtain the explicit notation of the mean molecular weight: $`\mu =[(4\pi G^3M^2)^1(k_\mathrm{B}/m_\mathrm{p})^4(3/a)]^{1/4}\{\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}^{1/2}`$. Taking the allowed parameter region into consideration, this can be expressed as

$$\mu (M,\beta _\mathrm{c})=1.06\left(\frac{M_{\mathrm{}}}{M}\right)^{1/2}\{\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}^{1/2}1,$$
(12)
which corresponds to the LE3 solution of $`\mu _{\mathrm{LE3}}=4.25(M_{\mathrm{}}/M)^{1/2}(1\beta _\mathrm{c})^{1/4}/\beta _\mathrm{c}`$. Furthermore, using equations (8) and (10), we obtain $`\rho _\mathrm{c}=\{(4\pi GR^2)^3[k_\mathrm{B}/(\mu m_\mathrm{p})]^4(3/a)\}^{1/2}\xi _1^3(\beta _\mathrm{c})`$. Making use of expression (12), this can be written as

$$\rho _\mathrm{c}(M,R,\beta _\mathrm{c})=4.73\times 10^{14}$$

$$\times \left(\frac{M}{M_{\mathrm{}}}\right)\left(\frac{100\mathrm{au}}{R}\right)^3\frac{\xi _1^3(\beta _\mathrm{c})}{\{\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}}\mathrm{g}\mathrm{cm}^3.$$
(13)
Thus, for given $`M`$ and $`\mu `$, the fraction of gas pressure $`\beta _\mathrm{c}`$ is identified by equation (12). In addition, for given $`R`$, the central mass density is determined by equation (13), in contrast to the LE3 solution of $`\rho _{\mathrm{c}\mathrm{LE3}}=7.68\times 10^{12}(M/M_{\mathrm{}})(100\mathrm{au}/R)^3\mathrm{g}\mathrm{cm}^3`$, which is independent on $`\beta _\mathrm{c}`$. These dependencies are summarized in Table 1.

In Fig. 1, we plot equations (12) and (13) as a function of $`(p_\mathrm{r}/P)_\mathrm{c}=(1\beta )_\mathrm{c}`$. It is found that, for $`M=M_{\mathrm{}}`$ as an example, the inequality of equation (12) requires the fraction of radiation pressure at the centre to take its parameter region of $`(1\beta )_\mathrm{c}6.7\times 10^3`$. 333For plasmas, replace with $`4.4\times 10^4(1\beta )_\mathrm{c}7.4\times 10^2`$, corresponding to the allowed parameter region of $`1/2\mu 2`$ (see footnote 1). If we specify the mean molecular weight, then the fraction can be fixed numerically. For $`\mu =2`$ (hydrogen molecule), we obtain $`(1\beta )_\mathrm{c}=7.4\times 10^2`$. In such a gas dominant regime of $`(1\beta )_\mathrm{c}1`$, the dimensionless quantity $`\rho _\mathrm{c}(R^3/M)`$ does not largely depend on $`\beta _\mathrm{c}`$, as seen in the figure. For $`R=100\mathrm{au}`$, the numerical solution reveals the central density of $`\rho _\mathrm{c}=8.6\times 10^{13}\mathrm{g}\mathrm{cm}^3`$, which is about $`90`$ per cent lower than that of $`\rho _{\mathrm{c}\mathrm{LE3}}=7.68\times 10^{12}\mathrm{g}\mathrm{cm}^3`$. For massive objects of $`MM_{\mathrm{}}`$, the radiation dominant regime appears. In the limit of $`\beta _c0`$, both $`\rho _\mathrm{c}`$ and $`\mu `$ are asymptotically close to $`\rho _{\mathrm{c}\mathrm{LE3}}`$ and $`\mu _{\mathrm{LE3}}`$, respectively. This characteristic is a result of the behavior of $`\mathrm{\Gamma }_14/3`$ for $`\beta 0`$ in equation (4).

In Fig. 2, for $`0.01(1\beta )_\mathrm{c}0.9`$, we show the radial profiles of the normalized mass density $`\rho (r)/\rho _\mathrm{c}`$, enclosed mass $`m(r)/M`$, total pressure $`P(r)/P_\mathrm{c}`$, as well as the adiabatic exponent $`\mathrm{\Gamma }_1(r)`$, comparing the LE solutions of index $`n=3`$ and $`3/2`$. Here the radius $`r=\alpha \xi `$ is normalized by $`R`$. It is shown in Fig. 2(a) that the density profile gradually deviates from that of the LE3 solution as $`(1\beta )_\mathrm{c}`$ decreases. Concerning the case that the parameters of $`M`$ and $`R`$ are fixed, in the lower $`(1\beta )_\mathrm{c}`$ the central density decreases, because the density varies sufficiently slowly from the centre to the envelope. In the limit of $`\beta _\mathrm{c}1`$, the present solution shows an exact agreement with that of the LE equation with a polytrope index of $`n=3/2`$. The radial profiles of the enclosed mass are shown in Fig. 2(b). For higher $`(1\beta )_\mathrm{c}`$ including the LE3, the mass is more condensed in the central region. For example, the radius confining the half-mass indicates $`R_{\mathrm{M}/2\mathrm{LE3}}=0.325R`$. As $`(1\beta )_c`$ decreases, this radius gradually shifts outward, up to $`R_{\mathrm{M}/2\mathrm{LE1}.5}=0.521R`$. It is confirmed that the gaseous polytropic effects are likely to flatten the density profile significantly.

In Fig. 2(c), we show the radial profiles of $`\mathrm{\Gamma }_1[\beta (\xi ),\gamma =5/3]`$, defined by equation (4). As the argument $`\beta (\xi )`$ monotonically increases outwards, the adiabatic exponent also monotonically increases, yielding the self-consistent variable polytrope. This is one of the most important results in the present paper. For all $`\beta _\mathrm{c}`$ being considered, the exponent varies quite slowly near the central region, while for $`(1\beta )_\mathrm{c}0.5`$, the variation around the envelope tends to be very steep.

In Fig. 2(d), we show the radial profiles of total pressure of $`P(\xi )=[GM^2/(4\pi R^4)]\{\xi _1^2(\beta _\mathrm{c})/\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}^2\{[1\beta (\xi )]/\beta ^4(\xi )\}^{1/3}\stackrel{~}{\rho }^{4/3}(\xi )`$, normalized by the central pressure $`P_c`$. It is noteworthy that, for the gas dominant regime of $`(1\beta )_\mathrm{c}1`$, the pressure formula exhibits the density dependence approximated by $`P(\xi )\stackrel{~}{\rho }^{5/3}(\xi )`$ with $`\mathrm{\Gamma }_1\gamma `$, fully consistent with Fig. 2(c). At $`\xi =0`$, the pressure takes the peak value, to give

$$P_\mathrm{c}=4.19\times 10^3\left(\frac{M}{M_{\mathrm{}}}\right)^2\left(\frac{100\mathrm{au}}{R}\right)^4$$

$$\times \frac{\xi _1^4(\beta _\mathrm{c})}{\stackrel{~}{m}^2[\xi _1(\beta _\mathrm{c})]}\left[\frac{(1\beta _\mathrm{c})}{\beta _\mathrm{c}^4}\right]^{1/3}\mathrm{ergs}\mathrm{cm}^3.$$
(14)
Notice that, in contrast with the LE3 solution of $`P_{\mathrm{c}\mathrm{LE3}}=0.582(M/M_{\mathrm{}})^2(100\mathrm{au}/R)^4\mathrm{ergs}\mathrm{cm}^3`$, equation (14) does depend upon $`\beta _\mathrm{c}`$, but weakly for $`(1\beta )_\mathrm{c}1`$. For $`M=M_{\mathrm{}}`$, $`\mu =2`$, and $`R=100\mathrm{au}`$, we obtain the central pressure of $`P_\mathrm{c}=4.6\times 10^2\mathrm{ergs}\mathrm{cm}^3`$, which is an order of magnitude smaller than $`P_{\mathrm{c}\mathrm{LE3}}`$.

Invoking the EOS, the temperature profile can be described as $`T(\xi )=[3GM^2/(4\pi aR^4)]^{1/4}\xi _1(\beta _\mathrm{c})`$ $`\times \{\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}^{1/2}\{[1\beta (\xi )]/\beta (\xi )\}^{1/3}\stackrel{~}{\rho }^{1/3}(\xi )`$ (not shown in figure). At $`\xi =0`$, it takes the peak value of

$$T_\mathrm{c}=1.14\times 10^3\left(\frac{M}{M_{\mathrm{}}}\right)^{1/2}\left(\frac{100\mathrm{au}}{R}\right)$$

$$\times \frac{\xi _1(\beta _\mathrm{c})}{\{\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}^{1/2}}\left[\frac{(1\beta _\mathrm{c})}{\beta _\mathrm{c}}\right]^{1/3}\mathrm{K},$$
(15)
whereas the LE3 solution reads $`T_{\mathrm{c}\mathrm{LE3}}=3.90\times 10^3(M/M_{\mathrm{}})^{1/2}(100\mathrm{au}/R)(1\beta _\mathrm{c})^{1/4}\mathrm{K}`$; both having the dependence of $`\beta _\mathrm{c}`$. For $`M=M_{\mathrm{}}`$, $`\mu =2`$, and $`R=100\mathrm{au}`$, the central temperature of equation (15) is found to be $`T_\mathrm{c}=1.1\times 10^3\mathrm{K}`$, which is again lower than $`T_{\mathrm{c}\mathrm{LE3}}=2.0\times 10^3\mathrm{K}`$. It is found that the decrease of the temperature is relatively small, within a factor of 2 for each $`\beta _\mathrm{c}`$.

### 3.2 Their application to the molecular cloud condensation: the critical radius of the ’first core’

In the context of the study of molecular cloud condensation in the ISM, the present quasi-stationary model is now available to provide an insight into the complicated dynamics (Penston, 1969; Larson, 1969; Shu, 1977; Saigo & Hanawa, 1998), in particular, the formation of the first core (Masunaga et al. 1998). For such an application, the effects of magnetic fields and turbulence might be taken into account, and the additional pressure $`\delta P`$ might be effectively included in equation (1) by replacing $`\mu `$ with $`\mu ^{}=\mu \beta _\mathrm{p}/(1+\beta _\mathrm{p})`$, where $`\beta _\mathrm{p}=p_\mathrm{g}/\delta P>0`$ (Bludman & Kennedy, 1996); e.g. for pure magnetic pressure, $`\beta _\mathrm{p}=(8\pi p_\mathrm{g})/B^2`$. Moreover, for $`T_\mathrm{c}2000\mathrm{K}`$, the dissociation of hydrogen molecules triggers gravitational contraction of the quasi-stationary first core. By invoking the scaling of $`T_\mathrm{c}1/\mu ^{}`$ and equation (15), therefore, we find the relation between the mass and radius of the first core as follows:

$$R13\left(\frac{1+\beta _\mathrm{p}}{\beta _\mathrm{p}}\right)\left(\frac{M}{0.05M_{\mathrm{}}}\right)^{1/2}$$

$$\times \frac{\xi _1(\beta _\mathrm{c})}{\{\stackrel{~}{m}[\xi _1(\beta _\mathrm{c})]\}^{1/2}}\left[\frac{(1\beta _\mathrm{c})}{\beta _\mathrm{c}}\right]^{1/3}\mathrm{au},$$
(16)
for $`\mu =2`$. Note that the constraint of equation (12) gives the dependence of $`\beta _\mathrm{c}=\beta _\mathrm{c}(M,\mu =2)`$ in equation (16). On the other hand, in the LE3 we obtain the mass-radius (M$``$R) relation of $`R4.6(M/0.05M_{\mathrm{}})\beta _\mathrm{c}\mathrm{au}`$, where $`\beta _\mathrm{c}=\beta _\mathrm{c}(M,\mu _{\mathrm{LE3}}=2)`$.

In Fig. 3 for $`\beta _\mathrm{p}1`$, we plot the numerical solution of equation (16), that is, the allowable parameter region in the M$``$R plane. For convenience sake, we also plot the fractional radiation pressure at the centre as a function of the mass of the core: $`(p_\mathrm{r}/P)_\mathrm{c}=1\beta _\mathrm{c}(M,\mu =2)`$. For $`MM_{\mathrm{}}`$, corresponding to $`\beta _\mathrm{c}1`$ for fixed $`\mu `$ (compare Table 1), the minimum core radius of equation (16) can be well approximated by that from the LE3: $`R_{\mathrm{min}}1.9\times 10^3[(1+\beta _\mathrm{p})/\beta _\mathrm{p}](M/100M_{\mathrm{}})^{1/2}\mathrm{au}`$. On the other hand, for $`MM_{\mathrm{}}`$, we numerically find that the minimum possible radius of equation (16) scales as

$$R_{\mathrm{min}}3.0\left(\frac{1+\beta _\mathrm{p}}{\beta _\mathrm{p}}\right)\left(\frac{M}{0.05M_{\mathrm{}}}\right)\mathrm{au},$$
(17)
and $`(1\beta )_\mathrm{c}2.8\times 10^4[\beta _\mathrm{p}/(1+\beta _\mathrm{p})]^4(M/0.05M_{\mathrm{}})^2`$. The estimation of (17) reasonably supports the results of radiation hydrodynamic calculations by Masunaga et al. (1998). They show that the mass and radius of the first core are $`M0.05M_{\mathrm{}}`$ and $`R5\mathrm{au}>R_{\mathrm{min}}`$, respectively, and the results do not largely depend on the initial conditions of wide parameter ranges. If $`RR_{\mathrm{min}}`$, the core tends to be gravitationally unstable and evolves along the dynamically-contracting track, to self-organize the ’second core’ as a protostar (Masunaga & Inutsuka, 2000).

## 4 Concluding remarks

In conclusion, we have developed a generic scheme to construct an exact polytropic sphere of self-gravitating gaseous and radiative medium. The particular results derived from the newly modified Poisson equation for hydrostatic equilibria are:

1. the numerical solutions show that for all cases, the central density and pressure are lower than those from the LE function with the index of $`n=3`$;
2. the adiabatic exponent monotonically increases radially outwards; and
3. in the special case that the central polytrope is gaseous closely, the present solution reproduces the properties of the LE function with $`n=3/2`$.

Within this framework, the whole mass of the system is connected with the central density, temperature, and the mean molecular weight.

For an application to modelling the molecular cloud condensation in the ISM, we have newly found the scaling law of the critical radius of the first core. The preliminary result is in consistent with that of the radiation hydrodynamic simulations. We expect that the major consequence can be also referred to, for example, mutatis mutandis, the study of the passive phase of protoplanetary discs (Honda & Nakagawa, 1999), stellar modelling (Basu, Pinsonneault & Bahcall, 2000), and so on.

## Acknowledgments

We acknowledge the Information and Educational Center, Kinki University TC for their hospitality. YSH is grateful to Yoshitsugu Nakagawa for a useful discussion.